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ABSTRACT 


This work concerns the partial structure factors of 

classical simple liquid mixtures near phase separation* The 

theory is developed for particles interacting through pair 

potentials, and is thus appropriate both to insulating fluids, 

and also to metallic systems if these may be described by an 

effective ion-ion pair interaction* The motivation arose 

from consideration of metallic liquid mixtures, in which 

resistive anomalies have been observed near phase separation. 

The pair potential is written as v.. = v . . + v. where v. . 

xj 13 3 .J 

is the pair interaction appropriate to a reference fluid and 

is chosen so that v}. may be treated as a perturbation. We 

J 

study how to correct a mean field theory appropriate to such a 

potential for the effects of correlated motions in the 

reference fluid. The work is cast in terms of functions f.. 

lj 

which are closely related to the direct correlation functions 
of Ornstein and Zernike. Exact equations for the f^ are 
derived by a method, originally developed for the study of 
the quantum. electron gas, which treats the densities of the 



mixture as basic variables in a linear response problem. We 
obtain approximate solutions to these equations which, however, 

I 

are exact to first order in . in the long wavelength limit, 
and which explicitly include the effects of reference system 
correlations. These solutions are then used to calculate the 
long wavelength form of the structure factors of metallic 
alloys, where we select for vY? f the potentials appropriate to 
a mixture of hard spheres. We seek to observe the singular- 
ities, at k » 0, associated with phase separation (in the 
critical region) and the long wavelength behavior which 
accompanies these singularities. The results are qualita- 
tively in accord with our physical expectations. Quantitative 
agreement with experiment seems to turn on the selection of 
the hard core reference potential in terms of the metallic 
effective pair potential, a task for which a successful 
systematic procedure has yet to be found. It is ^suggested 
that the present effective pair potentials are perhaps not 
properly used to calculate the metallic structure factors at 
long wavelength. Suggestions are made for application of 
these results to the thermodynamic and structural properties 
of insulating fluids. In the case of metallic systems, a 
qualitative explanation of the resistive anomaly is proposed, 
and suggestions made for a quantitative test of the hypothesis. 





I. Introduction and Statement of the Problem 

The work reported here arose from consideration of a 
single experimental result, the resistance anomaly observed 
in some binary liquid metal alloys as the phase separation 
temperature is approached. The importance of this result is 
that the theoretical difficulty it presents has its origin 
in a central problem of the theory of classical liquids, 
namely, the calculation of liquid structure from the 
interactions among the particles of the liquid. In this 
section, we shall describe these experiments and develop 
our analysis of these experiments to the point where a clear 
statement of the problem and program of this thesis can be 
made. 


I-A. Phase Separation and the Experiment of Schurmann 
and Parks 

To set the stage, we must begin with a recital of the 
basic facts of phase separation in binary systems. Consider 
a binary system composed of particles of type 1 and Nj 
particles of type 2. The thermodynamic state of this 

system may be considered a function of three variables? 

<*■ 

pressure p, temperature T, and concentration x = ^/(N^ + N^) • 
The phase diagram relevant to this work is that formed, at 
constant p, in the x-T plane. Fig. 1 represents schematically 
a portion of an x-T diagram of the type important to our work. 
It represents, in fact, an abstraction from the x-T diagram 

-i -b 



Concentration 


Figure 1 

Portion of an x-T phase diagram of a binary system 
which exhibits phase separation. 
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compiled by Schurmann (1971) for the Li-Na system. 

The phase diagram specifies, for each point (x,T) , the 
stable configuration of the corresponding mixture. In Fig. 1, 
points from the small triangular region at the lower left 
correspond to a solid phase. The phase boundary from D to 
E separates the remainder of the region shown into two 
liquid phases. For a point above the boundary, say the point 
A specified by ( X A * T A ) t the stable configuration is a single, 
uniformly mixed, alloy of concentration x A « For a point 
below the boundary, say the point specified by (x^ ,T^) , the 
stable configuration is a separated phase. In this phase, 
two distinct alloys are present simultaneously. (Under the 
influence of gravitation, these alloys will occupy two 
separate regions of the container, divided by a meniscus.) 
These alloys are described by concentrations x^ < x fi and 
x 2 > x b' w ^ ere x i x 2 are determined by T B alone as the 
concentrations at which the phase boundary cuts the line 
T = Tg. This construction is illustrated in Fig. 1. The 
relative amounts of these two alloys are determined from 
and x 0 , the given x_, and the requirement that all particles 
be accounted for. 

At points above the phase boundary, where a single 
uniform alloy can be maintained at equilibrium, the components 
are said to be miscible. The boundary defines a critical 
temperature T^ above which the components are miscible in all 
proportions. The point C at which the line T = T c intersects 
the phase boundary is the critical point, and defines a 
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critical concentration x c . T c and x Q are indicated in Fig. 1. 
The transition of an alloy from the uniform to the separated 
phase is known as phase separation. 

The experiments which motivated this work are concerned 
with what happens in the stable uniform alloy at temperatures 
just above phase separation. Schiirmann and Parks (1971) have 
measured the electrical resistance as a function of 
temperature at various concentrations in two metallic binary 
systems, Li-Na and Ga-Hg. As we shall outline shortly, the 
electrical resistance provides an integrated measure of 
fluctuation effects. The purpose of these experiments was 
to determine if changes in the fluctuation spectrum would 
cause the resistance to show any precursive behavior, that 
is, to respond in an observable way as phase separation is 
approached. Such a response might be considered analogous 
to the critical opalescence of the critical liquid-gas system. 

No precursive behavior was observed at any concentration 
in the Ga-Hg system. In the Li-Na system, some degree of 
precursive behavior appeared at most concentrations, and the 
effect became more pronounced as the concentration varied 
toward the critical value. The effect takes the following 
form. If the temperature is decreased, starting at a point 
well removed from phase separation, the resistance decreases 
at first in a nearly linear fashion. As the temperature of 
separation is approached, however, the resistance begins to 
decrease more rapidly, so that a plot of resistance vs. T 
develops a pronounced curvature, and the resistance attains 
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a value at separation which lies below the linear 
extrapolation of the high temperature results. 

I-B. The Resistance Anomaly and Liquid Structure 

To understand the significance of these experiments, we 
consider first the theory of conduction in liquid metal 
alloys. The simplest successful theory was proposed for the 
case of a pure fluid by Ziman (1961) . The theory was 
generalized to the case of binary alloys by Ziman and Faber 
(1965) . The development of the theory has been aided by the 
introduction of adequate models for the electron-ion 
pseudopotential and the liquid structure. Calculations for 
pure simple liquid metals (Ashcroft and Lekner, 1966) and for 
liquid binary alloys of simple metals (Ashcroft and Langreth, 
1967, A, B, to be referred to as Al-I and Al-II respectively) 
have achieved reasonable agreement with experiment. Although 
this success remains something of a mystery, we base our 
analysis on the Ziman theory. The theory, as we shall see, 
encompasses, in a reasonable way, the possibility of an 
anomaly near phase separation. 

In the formulation of Al-II, the result of the Ziman and 
Faber theory for the resistivity of a binary alloy is written 

p elec = Z * f dY y3{2/jrri ~ x) s i2 ( y )v i ( y )v 2 (y) 

e Jo 

+ xS 1]L <y) (Vj^ (y) ) 2 + (l-x)S 22 (y) (v 2 (y) ) 2 }. (1-1) 

For this discussion, we need not give careful definitions to 
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most of the pieces of this expression. The liquid metal is 

viewed as a system of rigid ion cores, of types one and two, 

moving about in a sea of nearly free conduction electrons. 

Defined relative to this model, kp is the usual Fermi wave 

vector, Z* is a (concentration dependent) effective valence, 

x is the concentration variable we have already defined, y is 

a wave number variable measured in units of 2k„, and the 

F 

functions v^(y) are dimensionless form factors for the 
screened electron-ion interactions. The key functions of 
our work are the functions (y) , the partial structure 
factors, for they are the only functions in (1-1) which 
contain detailed information on the ionic positions. They 
are defined by 

S ijQi) = (N i Nj)" Js <p i (k) 3 j (-k) > - (N^N..)* 5 6 k , 0 . (1-2) 

Here, and ^ are the particle numbers defined above, < > 
denotes a thermal average, and the operators ^ (k) are 
defined by 

ik * r 1 

e ± (k) = X e “ , (1-3) 

m 

i tVi 

where r^ denotes the position of the m particle of type i. 

The operators £^(k) have the property that 

P ± (k) = <3 i ‘(k) > / (1-4) 

th 

where (k) is the k Fourier coefficient of the equilibrium 
density of particles of type i. From (1-2) , we see that the 
<k) represent the static density fluctuation spectrum of 
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the liquid. It is reasonable to expect the density 
fluctuations to reflect in some way the nearness of the 
phase separation instability. This idea provided the 
original motivation for the resistivity measurements, and 
lies at the heart of our work. 

It will be useful in this work to present here a few more 
relations involving the structure factors. A more detailed 
discussion of such matters is reserved for an appendix (A) 
on the classical n-body distribution functions. 

Many of the general points we need to make can be made 
in reference to pure fluids, avoiding the clutter of notation 
introduced by consideration of binary systems. For a pure 
fluid, we can define a single static structure factor S (k) by 

ik •(r-r) 

S(k) =| <E e > - N <5 k Q , (1-5) 

nm 

where N is the number of particles and r^ denotes the position 
of the n particle. This function is related by Fourier 
transform to the radial distribution function g(r): 


g(r) 



e (S (k) - 1) 


( 1 - 6 ) 


where V is the volume of the system. The function g (r) is 
proportional to the conditional probability that a particle 
will be found at r, given a particle at the origin. Finally, 
the functions S (k) and g(r) are related to the two body 
distribution function .p 2 (r,r')- This function has the property 


*The integral in this equation, incidently, defines the 
Fourier transform convention which we adopt throughout this 
work. 
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that p 2 (r, r ' ) dr dr’ is the probability of finding simulta- 
neously one particle in the volume element dr about r and 
another in the element dr ‘ about r ’ . In the translationally 
invariant system, P 2 (*>£.' ) - p 2 (£-£.' ) and 

g(r) = (V/N ) 2 p 2 (r) . (1-7) 


The generalization of (1-6) to binary systems is 

r dk 


g. . (r) - 1 + 
^13 — N 


e x — *— (S . . (k) - 6 . .) , 

(2n) 3 13 ,3 


( 1 - 8 ) 


where g. . (r) now represents, suitably normalized, the 
probability of finding a particle of type j at a distance 
r from a given particle of type i. Clearly, g^j (r) = 9jj_(£)» 
The generalization of (1-7) is straightforward. 


We can collect what has been said so far into a first 
statement of the problem of this work. That problem is to 
calculate the partial static structure factors of the alloy 
with sufficient precision that the onset of phase separation 
can be observed and the resistance anomaly explained. 

I-C. Expectations 

The development of our approach has been guided to a 
considerable extent by our expectations for the form of the 
outcome. 

The first, and most important, of these expectations 
concerns the limit as k-> 0 , and takes the following form. 

At all points (x,T) , we can in principle calculate the 
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structure factors of the uniform phase. (For this purpose/ 
points below the phase boundary are exceptional only in that 
the appropriate Gibbs energy is higher than that of some 
separated phase.) A thermodynamic analysis (outlined in 
appendix B) shows that such calculations will yield, at each 
concentration, a temperature at which the partial structure 
factors will diverge as k->0. When plotted vs. concentration, 
these temperatures will form a curve suggesting, but not 
strictly reproducing, the phase boundary. This curve lies in 
general below the phase boundary, rising to join it at only 
one point, the critical point. We shall seek to calculate 
this line of hypothetical singularities at k - 0, and identify 
its highest point as the critical point. 

A second guide in this work has been a conjecture that 

it is the small k behavior associated with this k = 0 

singularity which gives rise to the resistance anomaly. As 

this conjecture is at variance with that of Schurmann and 

Parks (1971) , we should outline the manner in which the small 

* 

k behavior might account for the observed effects. First, in 
the resistance data for the Li-Na system, some degree of 


*This conjecture is also at variance with that of Fisher and 
Langer (1968) , who suggest that the observed resistance 
anomalies at magnetic critical points must not be due to long 
range correlations, because of the finite mean free path of 
the conduction electrons. Their point, however, is that the 
electron mean free path is limited by scattering other than 
spin scattering (e.g. phonons, impurities), while for the 
liquid metal alloy, we consider that the only scattering 
present is the scattering from density fluctuations treated 
in the Ziman formula. The electrons may thus be scattered 
by even long wavelength fluctuations. 
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an omalo us behavior is observed at most concentrations, but 
the effect clearly becomes more pronounced as the concen- 
tration is varied toward the critical value. This observation 
seems to reflect the relation, described above, between the 
phase boundary and the line of singularities. 

Secondly, we must ask why the effect is observed in 
Li-Na but not in Ga-Hg. Consider the form of (1-1). The 
upper limit on the integral is 2k p , which depends on the 
effective valence of the alloy. Now the partial structure 
factors are characterized by strong first peaks (AL-II) , and 
the valence dependence of the integration limit is such that 
these peaks lie within the range of integration for Ga-Hg, 
but just outside it for Li-Na. Thus, when we observe that 
large k effects are heavily weighted by the factor y 3 , it is 
reasonable that a small k effect will be lost in Ga- Hg, but 
(barely) observable in Li-Na. We note in support of these 
ideas that the effect in Li-Na is very small. Even at the 
critical concentration the resistance at separation lies less 
than 1% below the linear extrapolation of the high temperature 
results. 

Finally, we must ask why the resistance is depressed by 
this effect, when the structure factors of (1-1) are diverging. 
When Bhatia and Thornton (1970) discussed these singularities, 
they in fact suggested that the resistance should go up. In 
this matter, we can only note that the complexities of the 
alloy (e.g. the factors and S 22 diverge positively while 

S ^2 diverges negatively) are such that a clear prediction may 
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not be possible. 


I-D. Structure of Metallic Liquids - Previous Results 

To arrive at a final statement of the problem, we need 
to review the assumptions and principle results of the theory 
of liquid metal structure, as it existed at the inception of 
this work. 

At the level appropriate here, the theory of liquid 

metals rests on several standard assumptions. The fundamental 

assumption is the validity of the adiabatic approximation, 

which asserts that, on the time scale of ionic motions, the 

relatively light and mobile conduction electrons adjust 

instantly to any change in ionic configuration. This means 

that, for purposes of determining ionic motion, the electronic 

configuration may be considered to be completely specified by 

only the volume V and the ionic positions, which we denote by 

the general variable R. Then the total energy may also be 

considered a function of only V and R, and denoted by 

E metal^ V,R ^* A second assumption (actually implicit in this 

definition of E . . (V,R)) is that the ions exhibit no 

metal 

internal structure. The final assumption is that the ions 
form a classical system so that their motion may be determined 
by classical dynamics. 

With these assumptions, the problem of ionic structure 
factors and motion is cast as the problem of a classical 
liquid, with, however, a complicated energy function. In 
theories of classical liquids, the assumption is almost 
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uni vers ally. made that the energy function depends only on the 
particle positions R, and is expressible as a sum of pair 
interactions. That is, it is assumed that the energy may be 
written as 

U(R) = | C v(r. - r.) . (1-9) 

Ashcroft and Langreth (1967C, AL-II) have proposed an 
approximate form for E metal < v / R ) in which all dependence on 
R can be accounted for in a sum over an effective, volume 
dependent, pair interaction between ions. (The total energy 
of the liquid, however, is not recovered by a sum of the form 
(i“9 ) t since the sum excludes some volume dependent terms.) 
This effective pair potential will be denoted by v AL . These 
results suggested that classical liquid theory might be 
applied to the ions of a simple liquid metal at constant 
volume by using v AL as the pair interaction. This idea has 
opened the way to important progress in understanding liquid 
metal structure. Though we shall mount a selective challenge 
to this idea in section V, it formed a crucial starting 
assumption of the present work. 

The important strides in understanding liquid metal 
structure have come with the application of the hard sphere 
model. When the potential v AL is viewed in real space (AL-II) , 
it is characterized by a short range, harshly repulsive core, 
and a weak, long range, generally attractive tail. The 
simplest non-trivial model of such a potential is that acting 
between hard spheres of diameter c: 



The structure factors appropriate to the hard sphere model 
fluid are reasonably well known from machine calculations 
(Alder and Wainwright, 1958; Alder and Hoover, 1968; Wood, 

1968; Rowlinson, 1968) and in the analytic Percus-Yevik (PY) 
approximation (Percus and Yevik, 1958; Wertheim, 1963; Thiele, 
1963; Lebowitz, 1964; Ashcroft and Lekner, 1966; AL-I ). 

Even before the form of v^ was suggested, Ashcroft and 

Lekner had applied the hard sphere model to pure simple liquid 

metals with considerable success. They showed that the hard 

sphere structure factors of the PY approximation could be 

adjusted by selection of a to give an excellent fit to the 

experimental structure factors of these systems in the region 

of the first peak and, to a lesser extent, down to the 

smallest wave vectors for which experimental data were 

available. Using these hard sphere structure factors in the 

Ziman formula, they calculated the electrical resistivity and 

achieved reasonable agreement with experiment. They also 

noted that the hard sphere diameter a used to achieve the 

structural fit corresponded, near the solidification point of 

* 

each element, to a packing fraction of about 45%. This 
compares favorably with the value of 49% at which machine 
calculations place the crystallization of the hard sphere 
liquid (Hoover and Ree, 1968). 

*The packing fraction is the ratio of the volume occupied by 
the hard spheres to the total volume of the system (eq. (5-9)). 
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AL 

When the potential v was introduced, Ashcroft and 
Langreth showed that it bore a simple physical relation to the 
hard core diameters of the Ashcroft-Lekner work. They went on 
to generalize this relation to the case of a binary alloy, and 
use the structure factors of the hard sphere mixture (in the 
PY approximation) in a successful calculation of alloy 
resistivities (restricted, of course, to alloys which do not 
separate) . 

The hard sphere success, however, does not appear to 
include an accounting for phase separation. Phase separation 
has not been observed in hard sphere mixtures either by 
machine calculation (Alder, 1964 (M.D.) and Rotenberg, 1965 
(M. C . } ) or by analytic study from the PY approximation 
(Lebowitz and Rowlinson, 1964). Although this is apparently 
still an open question, the negative result is not unexpected. 
The hard sphere potential defines no energy scale, but only 
a scale of length. Then a phase transition in a hard sphere 
fluid can be driven only by geometry. When we consider phase 
separation in mixtures, we note that both the uniform and 
separated phases are characterized by disordered alloys, so 
that we do not expect geometry to play the dominant role in 
the transition. (We shall, in any case, treat the important 
hard sphere systems of this work in the PY approximation, 
within which these systems are known to be stable.) The plan 
of our work will then be to go beyond the hard sphere model in 
order to understand phase separation. 
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I-E. Statement of the Problem and Outline of the Thesis 

What is needed is a theory of structure which preserves 
the successes of the hard sphere model, but goes beyond this 
model in order to account for phase separation. The natural 
suggestion is that we try to calculate the liquid structure by 
a perturbation theory which begins in zero order with the hard 
sphere liquid. 

The actual program set for this work results from 
combining this suggestion with one further hypothesis. That 
hypothesis addresses the following question. The perturbation 
theory we envision will find its most natural formulation for 
potentials of the form 

v = v HC + v 1 , (1-11) 

where v 1 is the perturbation. But the potential v ,AL cannot be 

cast in this form, since all potentials of this form have a 

rigorous hard core, while v AL has a core with a steep but 

finite slope. It is our hypothesis that the "softness" of 

the core has little effect on the physics of phase separation, 

and thus, at least for a first attempt at this problem, we 

need not consider the more elaborate perturbation theory 

* 

required to treat this softness. Rather, we assume that the 
essential physics of phase separation is implicit in a 
potential of the form (1-11) , and develop our theory for such 
a potential. We call this potential the "hard core plus tail" 

*Such theories have been developed. See for instance Andersen 
et al. (1971). and Barker and Henderson (1967B) . 
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V AU V HC 



FIGURE 2 
AL HC 

Schematic showing v (r) and v . A hard core plus tail 

AL HC AL 

model of v will be v plus that part of v which lies 


outside the core. 
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potential. In Pig. 2 , we present a schematic showing a 
potential v^, and a hard core plus tail model of this 
potential. 

Based on these considerations, the program set for this 
work was to develop a suitable perturbation approach to the 
hard core plus tail liquid, and then to perform calculations 
for real liquid metal systems by constructing hard core plus 
tail models of the potentials v^. In section II, we discuss 
the general nature of the perturbation problem in liquids , and 
then present and discuss the mean field theories which 
provided the motivation for our approach. In sections III and 
IV, we develop our approach to the problem, first through a 
set of general equations, and then through a study of 
approximate solutions to these equations. The best solutions 
appear at long wavelengths. In section V, we turn to the 
application of these results to liquid metallic systems. We 
perform calculations for both pure fluids and binary systems. 
This work ultimately encounters a fundamental difficulty in 
the use of the pair potential v^ for calculations at long 
wavelength. We do not undertake to resolve this difficulty, 
but we do suggest that our calculations are perhaps best 
viewed as model calculations, which incorporate some, but not 
all, of the essential features of metallic systems. The 
calculated results are in accord with our physical expec- 
tations. Although a resistivity calculation is not presented, 
suggestions are made for the form of a calculation to test the 
hypotheses described above. Finally, in section VI, we review 
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and discuss what has been learned, and suggest further lines 
of work both for the perturbation theory itself and for its 
application to metallic systems. 



IX. General Discussion and Mean Field Theory 

In this and the next two sections, we shall develop our 
approach to the perturbation theory of the hard core plus tail 
mixture. In recent years, a good deal of effort has been 
devoted to the perturbation theory of simple liquids. (For a 
recent review, see Barker and Henderson, 1972). Some of this 
work, dealing with the softness of the core, has already been 
mentioned. The major lines developed in this work have been 
directed, initially at least, at the calculation of thermo- 
dynamic and structural properties of Lennar d- Jones fluids away 
from any critical points. In contrast, the thrust of our 
effort has been to develop an approach specifically designed 
for the problem of critical fluctuations in binary metallic 
liquid systems. The result has been that our work did not 
grow directly from any of the major lines of recent 
perturbation theory, but developed instead out of consider- 
ation of the simplest mean field theory of the hard core plus 
tail liquid. In this section, we shall first review briefly 
some of the main lines of liquid perturbation theory, for the 
insights they offer into the general nature of the problem, 
and then turn to a more detailed presentation and discussion 
of the mean field approach. In order not to clutter this 
section with the notation of mixtures, we shall make its 
general points in reference to pure fluids, except where 
explicitly noted. 
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II -A. General Discussion 

We seek a perturbation theory which begins with a 
"reference" fluid described by a pair interaction v re ^, and 
considers the properties of a "real" fluid described by a pair 
interaction, v, of the form 

ref 1 

v = v + v . (2-1) 

The plan is to calculate the properties of such a fluid as an 
expansion in v"* - about the properties of the reference fluid. 

We shall ultimately select for the reference potential the 
hard core potential of eq. (1-10). 

This program faces a major difficulty in principle, 
namely, that the ultimate reference system, the hard sphere 
liquid, is incompletely understood. The limits of our 
knowledge are best expressed in the language of the n-body 
distribution functions. (See Appendix A). Of these, only the 
one and two-body distribution functions are known with any 
accuracy for the hard sphere liquid. But straightforward 
expansions in of either the Helmholtz free energy (Zwanzig, 
1954) or structure factors (Coopersmith and Brout, 1963? also 
Brout, 1965) of the "real" fluid lead rapidly to terms 
requiring higher distribution functions of the reference 
fluid. This difficulty is of crucial importance in deter- 
mining the form ultimately taken by any of the perturbation 
theories, including our own. 

This point is usually made in another way. Although the 

V*0 f 

calculation of p ^ for the hard core liquid is a field in 
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itself, the fundamental assumption of all the perturbation 

rgf 

theories, again including our own, is that p 2 is a known 
function. 

The present activity in the perturbation theories of 
liquids begins with the work of Zwanzig (1954) who derived 
the expansion for the free energy of a fluid with the pair 
potential of the form (2-1) . The expansion takes the form 

F = F ref + jjdr dr’ v 1 (r-r 1 ) p 2 ® f (r ,r ' ) + ofv 1 ) 2 . (2-2) 

Here, F and F are the free energies of the "real" and 
reference fluids respectively, and the other elements have 
already been defined. All higher terms in this series involve 
reference distribution functions of order higher than two. In 
early work (e.g. Smith and Alder, 1959) , the series was 

3T0 f 

truncated after the first order term so that only p 2 was 
required. Though we start from a different viewpoint, an 
expansion similar to (2-2) eventually plays an important role 
in our work. 

In later work. Barker and Henderson (1967A) studied the 
second order term in (2-2) . This term contains reference 

2T6 f 3T©f 

three end four-body distribution functions, e 3 and p 4 , 
which, for these purposes, are essentially inaccessible to 
present analytic and "machine experiment" techniques. Viewed 
in one way, the work of Barker and Henderson shows that these 
complicated functions actually contain more information than 
is needed to calculate the second order term. They cast the 
required information into a form which is reasonably, 
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accessible to analytic approximation (1967A) and machine 
calculation (Barker and Henderson, 1972). The Barker- 
Henderson theory has found its chief application to the 
thermodynamic (and some structural) properties of Lennard- 
Jones fluids. (These are reviewed by Barker and Henderson, 
1972.) Though the approach and application are very different 
from ours, there are ultimately some points of contact between 
the Barker-Henderson work and our own. These will be pointed 
out in the concluding discussions of section VI. 

Turning from thermodynamics to consideration of the 
structure factors, we find that the direct expansion of S (k) 
in powers of v^ yields, even in first order, terms requiring 
hard sphere distribution functions of order greater than two. 
This expansion was presented by Brout and Coopersmith (1963) , 
who attempted to surmount this difficulty by approximating the 
higher order distribution functions with a superposition of 

ITS f 

the functions p 2 . Since their work, it has been shown that 
this problem has a solution in principle, which can be 
obtained by beginning the perturbation theory not about the 
hard core fluid, but about the ideal gas, for which, of 
course, all distribution functions are known. The result is 
a series in the full potential v. (The development of this 
series is presented in several places. For presentation, 
discussion and references, see Rushbrooke (1968).) If the 
potential v is then split into v and v , partial summations 
can be performed which eliminate v re ^ in favor of p~ ef , 

z 

yielding a series in which each term can be evaluated knowing 
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only and v 1 . This development forms the basis of recent 

work by Andersen and Chandler (1972) , to which the reader is 
referred for discussion and references. Though we shall not 
have occasion to use this series explicitly, we have in fact 
studied it in some detail, and its existence and form have 
influenced this work at some points. These points will be 
mentioned as they appear. The Andersen-Chandler work itself 
has some relation to our own, and we shall have occasion to 
mention it later on. 

As discussed at the outset, our work, despite its links 
to these various approaches, has developed primarily from 
consideration of the mean field theory of the hard core plus 
tail system, to which we now turn. 

II-B. The Mean Field Theory 

At the base of our work is the observation that the 
structure factor S(k) bears a simple relation to the static 
density response function of the liquid. Consider a static 
external field $ (r) applied to a pure liquid of mean density 
p 0 . The liquid will respond to this field by assuming a 
configuration in which the local density, given by p (r) , will 
not in general be uniform. If the density response 6p(r) = 
p(r) - p 0 is given in reciprocal space by 

6 p (k) = x(]0<K]0 + o(<)>) 2 , (2-3) 


then, for a classical system, 
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x(k) = -SPoS(k) . (2-4) 

(Here, 3 = l/k^T.) This result may be straightforwardly 
derived from the definition of S (k) and the statistical 
definition of the density thermal average. The function x (k) 
is the linear density response function of the liquid. Though 
the relation (2-4) is a simple proportionality, it is 
convenient to cast the development of this and the following 
two sections in terms of x rather than S, for that development 
leans heavily on the linear response interpretation of S. 

To display the mean field theories , and for later use in 
the formalism, it is useful to develop the following notation. 
We first give a name to the inverse of x by defining 


f (k) 


1 

X (k) * 


(2-5) 


The comparison of "real” and reference fluids then takes the 

following form. The "real" system, with potential v, is 

characterized by a response function x (k) t while the reference 

system, with potential v , is characterized by the response 
re f 

function x (k) ■ These functions in turn, through the form 
(2-5), define functions f (k) and f ref (k). Then, with the 
natural definition 


fl (k) - f(k) - f ref (k) , 
we can write 

ref ,, . 

X (*> 

1 - X ref (k) fhk) 


( 2 - 6 ) 


X(k) = 


(2-7) 
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If we substitute (2-7) into (2-3) , the result can be written 

6p(k) = x ref (k) 4 > eff ( k ) + o (*) 2 , ( 2 - 8 ) 

ef f 

where, in real space, the "effective potential" <J> takes the 
form 

* eff (r) = (r) + Jdr* f 1 (r-r ' ) Sp (r ' ) . (2-9) 

In words , these last equations cast the response of the real 
fluid as the response of the reference fluid to an effective 
potential which contains both the applied potential and the 
effects of the perturbing pair potential v\ 

The simplest mean field approximation is to let the 

reference system be an ideal gas (v =0) and take for <}> 

* 

the "Hartree" potential; i.e. 

X ref (k) = -8 Po (2-10) 


and 


<{> eff (r) = 4>(r) + Jdr' v (r-r ' ) «P (r ' ) 


Then, clearly f = v, and we have 


( 2 - 11 ) 


< 2 - 12 

These steps represent a classical realization of the familiar 
random phase approximation to the response of the quantum 


*The result (2-10) can be calculated directly, but it also 
follows quickly, via eqs. (1-6) and (2-4) , from the well known 
result that g(r) = 1 for an ideal gas. 
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electron gas. Stroud (1973) has already used the general- 
ization of this approach for a two component system to 
consider phase separation in binary liquid metal alloys. 

It is our purpose to consider the case in which v ref 
HC 

- v . A natural approximation (in the spirit of (2-11)) 
would be to let <£ e: ^ be a Hartree potential based only on v^: 

* eff <r) = <f> (r) + Jdr 1 v 1 (r-r ’ ) <$p (r ’ ) . (2-13) 


. 1 
In this form, f 


and we have 


XUO 


1 


x HC 


(k) 


HC ... 1 „ . * 

X (k)v (k) 


(2-14) 


Such a theory, however, suffers from a peculiar ambiguity. 

HC 

The problem is that, for v = v , the region r<a can never 
be sampled. Then the physics of the liquid must be indepen- 
dent of the form given to v (r) for r<a. This requirement is 
clearly not met by the form (2-14) . 

To perform, with any confidence, a modification of (2-14) 
which removes this difficulty , we must return to the beginning 
and develop a systematic theory. But that is not the purpose 
of this section. Our purpose here is to present the flavor 
of the mean field idea, and to describe some calculations 
which, by pointing out the faults and virtues of that idea, 
motivate the development of our formalism. For these purposes, 
the ambiguity of (2-14) is really only a technical point. We 
therefore proceed at this point simply by stating that if one 
considers the systematic perturbation series and the partial 
summation represented by (2-14) , it is seen that a simple 
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modi f i cation , which represents a similar partial summation, is 
to write 


'-MF 


X HC 00 

1 - x HC QO fi , <*> 


(2-15) 


where 


E - (r> = 


r<c 

r>a 


(2-16) 


In this form, the theory is independent of v A (r) for r<a. As 
Indicated in the notation, we shall refer to this as the mean 
field theory of the hard core plus tail liquid. 


II-C. Discussion of the Mean Field Theory 

Mean field theories represent a common first approach to 
the physics of a phase transition. In this problem as well, 
the mean field theory achieves a measure of success. In a 
calculation best detailed in section V, we applied the mean 
field approach to the Li-Na binary system. The program was 
to model the Li-Na system with a hard core plus tail mixture, 
and apply the generalization of (2-15) and (2-16) for this 
case. We sought to identify the locus in the x-T plane of 
points at which the mean field partial structure factors 
diverge as k approaches 0. Within a crude but reasonable 
model , the line of singularities can readily be calculated. 
The result is a curve which, in shape, symmetry, and position 
in the plane, strongly suggests the experimental Li-Na phase 
boundary. (See Fig. 3.) 

Despite this success, the mean field theory suffers from 
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a serious inadequacy when applied to such dense systems. The 
nature of the difficulty can be illustrated by considering an 
extreme case and then observing parallels between this case 
and the case of the dense fluid. Suppose we let the density 
increase until the hard core reference "fluid" forms a 
rigorously close packed lattice. Then the subsequent addition 
of any tail potential v^ will clearly have no effect on the 
structure, because the particles are unable to move in 

response to that potential. For any v^, we have in this 

HC HC 

limit x = X ♦ Now, the structure factor and hence x °f a 

close packed lattice of hard spheres is a series of 6-function 

spikes. The formula (2-15) is thus clearly inadequate. The 

origin of the difficulty is that the mean field formulation 

treats the response of each particle to the fields of the 

other particles as if those fields were part of the applied 

external field. In fact, the accessible responses to the two 

types of field are very different. In this limit of close 

packing, the particles will execute a strong collective 

response to an external field with the periodicity of the 

lattice (this is the meaning of the spikes in x) / but can 

execute no relative motion, and hence no response to fields 

fixed to the particles themselves. 

The essential features of this situation survive, in 
muted form, when we allow the density to relax to liquid 
densities. The particles can now execute some relative 
motion, but the motion is nevertheless severely limited by the 
hard core packing. In this sense, the success of the hard 
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sphere model at these densities (and away from critical 

points) is a reflection of the rigorous close packing result 
HC 

X = x • In reflection of another close packing result, the 
hard core structure factors at these densities are charac- 
terized by strong (but now finite) peaks, indicating strong 
collective responses to applied fields of certain wavelengths. 
Thus, when the effect of a tail potential is examined in the 
mean field theory, this theory erroneously predicts a large 
effect at the wavelengths of the hard sphere peaks. We 
observed this failure in a calculation of a hard core plus 
tail model of pure liquid liquid Na. For a model appropriate 
to temperatures just above solidification, the mean field 
theory failed catastrophically at the first peak. This 
failure completely destroyed the excellent agreement in this 
region between the structure of the simple hard core model 

and the experimental structure factor of Na. This calculation 

* 

is reported in section V. 


*An approach which appears to surmount this difficulty is the 
"optimized random phase approximation" of Chandler et al 
(1972) (see also Andersen and Chandler (1972)). This work 
replaces the mean field formulation (2-16) with the two 
equations 

f 1 (r) = v 1 (r) for r>a , 

and 

g (r) =0 for r<g 

leaving f ^ (r) for r<a to be determined. Here, g(r) is the 
radial distribution function, defined in (1-6), for the “real" 
fluid, so that the second equation represents an exact result 
for a hard core plus tail fluid. These equations turn out to 
be sufficient to determine f 1 (r) for all r, and hence x and 
S(k). The resulting structure factors (for a hard core plus 
tail model appropriate to a Lennard-Jones fluid) show quite 
reasonable behavior in the neighborhood of the first hard 
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Based on these considerations, an important thrust of our 
work was to discover how to treat these collective or 
"correlation" effects. Keeping in mind the limit of close 
packing, we sought a theory which would be adequate in this 
limit, and thus deal at least with reference system 
correlations. 

Though the mean field theory suffers from this difficulty, 
its success with the line of singularities at k = 0 has been 
an important influence on our work. First of all, it tends to 
support our hypothesis that the hard core plus tail potential 
contains the essential physics of phase separation. In fact, 
during the development effort, our faith in this hypothesis 
rested in large measure on a crude form of this calculation. 

A second influence is more important, and also more 

subtle. The point concerns the relation (2-7) between x» 
ref 1 

X , and f . This relation is, of course, just definition, 
and yet it represents a rearrangement of the perturbation 
series in the sense that x contains terms of all orders in f 1 . 
Thus, the mean field theory may be viewed equivalently either 
as an (incomplete) low order approximation to f^ or as a 
selection of terms of all orders in the series for x* Viewed 
in the light of this observation, the success of the mean 
field theory in locating the phase boundary suggests that, in 

(continued from previous page) sphere peak. The first formula, 
however, is quite arbitrary at this stage, and, as it repre- 
sents a severe restriction on the form of f 1 , considerable 
justification is required before the results can be understood. 
The systematic viewpoint and exact results developed in this 
thesis might possibly shed some light on its meaning. 
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some sense, whatever physics attaches to the function may 
be more simply related to the phase transition that the 
physics attached to x* Such considerations have implicitly 
motivated the shape of our effort, which is to study the 
physics of f\ and ultimately to seek a theory of phase 
separation by completing the low order approximation to f^. 

Our search into the physics of f 1 began with considera- 
tion of its role in the effective potential formulation 
defined by (2-8) and (2-9). By these definitions, the real 
fluid response to the applied field is cast as the reference 
system response to an effective field in which f^ plays a 
central role. We might ask, what is the correct effective 
potential? As it turns out, this question has been asked 
before, in the study of the dielectric response of the 
quantum electron gas. There, the reference system is a non- 
interacting Fermi gas, and the effective potential is usually 
formulated as 

<t eff = <b + d>„ + d> + <J> , (2-17) 

v Y Y H T ex Y corr. 

where <£„ is the Hartree potential of (2-11) , and <f> and 

£i GX 

d> are corrections for exchange and correlation. (See, 

corr. 

for example, Ballentine (1967).) A study of this work led us 
to study the work of Hohenberg and Kohn (1964) and Kohn and 
Sham (1965) , who, considering the non-uniform electron gas, 
presented a formally exact prescription for the potentials <f> 

wA 

and <b . Elements of their work, transcribed to the case 

corr • 
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of a hard core plus tail mixture, form the basis of our work. 
We now turn to the development of our approach. 



III. The Formal Development 
III-A. The Basic Equations 

In this subsection , we use the statistical mechanics of a 
classical fluid in the presence of external fields to derive 
the basic equations of our work. Because we shall need the 
generalization to binary alloys, we shall work from the outset 
with a multi-component mixture. 

We consider an m-component mixture containing parti- 
cles of type i , for i - l,2,...m. We shall evaluate thermal 
averages in a canonical ensemble at constant temperature T, 
volume V, and particle numbers N^. We consider the mixture in 
the presence of a set of external fields. The configurational 
free energy of the mixture may be written 


pU 

•4 


e ln n, in, 


, ! . . . N I J' 


dR e -* U(R) e , ( 3 - 1 ) 


1 *"2 m 

where the integral over dR is over all co-ordinates of the 
particles, U(R) denotes the potential energy of interaction 
among the particles, and 4 (R) represents the energy of inter- 
action between the particles and the external fields. We 
shall ultimately take U(R) to be a sum over pair potentials, 
though, for now, it can remain unspecified. For 4 (R) , we 
intend from the outset a sum over single particle potentials. 
That is , 


N, 

*<R) = £ if / 

i SL=1 1 * 


(3-2) 


where <Jk is a field which couples only to the particles of 
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type i, and denotes the position of the particle of 

type i. We shall ultimately consider the effects of varying 
both U and 4, and have thus included their specification , as 
super- and subscripts, in our notation for the free energy. 

We are interested in the linear response of the densities 
of this system to the external fields. The necessary termi- 
nology is defined in generalization of that describing single 
component fluids. We denote the local thermal average density 
Of particles of type i by p^ (r) . In the absence of any 
external fields, this density will be uniform, taking every- 
where the value p? = N^/V. When the external fields <j> are 
applied, the density need not be uniform. We denote the 
density in the presence of the fields by pf (r) , and define the 
density response by 

6pJ(r) = pj(r) - p? . (3-3) 

The most general linear relation between the 6pf and the cf>^ is 
of the form 

<$pt(k) - C X ±j QiHj (k) + o(<J>) 2 . (3-4) 

This relation defines the linear response functions x^j °f the 
mixture. A calculation analogous to that giving (2-4) yields 

Xij(k) = -e/p^pj S ±j (k) , (3-5) 

where the are the partial structure factors of (1-2). 

The following development parallels the work of Hohenberg 
and Kohn (1964) and Kohn and Sham (1965) , as generalized to 
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finite temperatures by Mermin (1965). As the argument is 
somewhat involved, it may be useful to give an overview before 
beginning in earnest. The physical point which provides the 
basis for the work concerns the relations between the den- 
sities (r) , the external fields (r) , and the free energy 
F^. We shall throughout the first part of the argument 
consider U(R) to be specified. Then it is clear that the 
potentials cJk are sufficient to determine the (r) and F^. 
This fact has been indicated above by placing the superscript 
$ on p^. Suppose instead we specify the densities p^ (r) , and 
denote by a set of external fields which gives rise to 
these p^(r). The important question is, to what extent are 

and hence F^ determined by the p . ? The answer can be 

1 

obtained from a variational principle and is, not surpris- 
ingly, that is determined to within constant terms, and 
hence, by (3-1) , F^ is also determined within constant terms. 
We can show further, however, that even the constants cancel 
from the combination 

G(p,U) = F U - E jar $ P (£) Pi (r) , (3-6) 

so that, as indicated in the notation, G is uniquely deter- 
mined by the densities p^ (r) (assuming, as stated above, that 
U(R) has been previously specified.) Then for given U, the 
function G(p,U) may be expanded (formally at least) in a 
functional Taylor series about some specified density 
functions p^ (r) . If the expansion is carried out about the 
particular set p? defined above, the coefficients of the 
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second order terms turn out to be just generalizations of the 
functions f introduced for pure fluids in section II. By 
writing G(p,U) for real and reference fluids, we thus generate 
equations for the functions which are generalizations of f"*". 
These form the "basic equations" promised in the title of this 
subsection. 

We shall, in this argument, consider at some points that 

the densities p (r) are the starting point, so that we 

i 

consider p^ and and at other points that the potentials <(k 

are the starting point so that we consider and pt. We 
shall endeavor to be clear, both in context and notation, 
which view we adopt at each point. At the outset, until we 
prove the uniqueness theorem, we must consider that we start 
with the potentials 

The argument begins with a variational principle. 

Consider the functional of P (R) , U(R) , and 0 (R) defined by 

= JdR P(R){U(R) + $(R) + In P(R) 


+ | C ln(N. 1) } , 

P ■; 1 


(3-7) 


where U and $ are defined above, and P (R) is some distribution 
function satisfying 

JdR P(R) = 1 . (3-8) 

If we select for P the particular function P^ defined by 

P*(R) - e ' eU(R) e -^ (R) /(J’dR e" eU(R) e ~^< R) ) , (3-9) 
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the functional ft reduces to the free energy of (3-1): 

Q(P*,U.,«) = . (3-10) 

A straightforward adaptation to a canonical ensemble of a 
grand canonical argument by Mermin (1965) yields the minimum 
principle 

ft(P,U,$) > n(P*,U,®) , (3-11) 

for any P ^ P^ which satisfies (3-8) . 

In using this minimum principle, we will limit our 
consideration to the set of P’s which may be generated by 
evaluating (3-9) for all possible fields 0. Since the special 
function P^ of the minimum principle is a member of this 
restricted set, the minimum principle will still apply within 
this restricted set. 

We can now prove the necessary uniqueness theorem. 
Consider two external fields $ and of the form (3-2), with 
single particle potentials <Jk (r) and (r) respectively. The 
uniqueness theorem asserts that if pt(r) = pf (r) for all i, 
then for each i, (r) - (r ) is a constant independent of 
r. To prove this, we suppose the opposite and deduce an 
absurdity. We suppose that pt (r) = pt (r) for all i, and that 
there exists at least one i for which ^(r) - <|>| (r) is not 
just a constant. The potentials 0 and $' will define, through 
(3-9) and (3-1) , distribution functions P^ and P^ and free 
energies and F^, . Under the hypothesis that 4 k (r) - <J>| (r) 
is not just a constant for at least one i, we have P^ ^ P^ , 
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so that the minimum principle applies. In particular 


fJ, - ,U,$' ) 


< = fl(P*,U,#) + JdR P(R) {*' (R)-$(R) } 


(3-12) 


or 

fJ. < F $ + E jdr pt(r){(j,: (r)-^i(r)} , (3-13) 

where, for the last step, we make use of (3-2) and the 
statistical definition of (r) given in appendix A. Now the 
argument to (3-13) can be repeated with primed and unprimed 
interchanged, yielding 

F $ < fJ, + E jdr pt' (r){* ± (r)-^(r) > . (3-14) 

If we now suppose that p?{r) = pt (r) for all i, addition of 
(3-13) and (3-14) yields immediately 0 < 0. This proves the 
•uniqueness theorem. 

With this theorem, we can now consider the meaning of 
taking the p^ (r) as our starting point. We take these to be 
any densities which may be induced in the system by some 
external fields 0. If we denote by 0 P any potentials which 
give rise to these densities, the uniqueness theorem asserts 
that $ p is determined within constants. If we then proceed to 
form a distribution function P p by inserting one of the 
potentials 4> p into (3-9) , we see that the undetermined con- 
stants in $ p cancel, so that the distribution function P p is 
determined uniquely by the p^ (r) . Then, by inserting P p into 
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(3-7) , we may define a unique functional of p, and U, 
denoted by F(p,$,U), as 


F(p,$,U) 5 ft(P p ,$,U) 

= JdR P P (R){U(R) + $(R) + ~ In P P (R) 


(3-15) 

+ In (N i ! ) } . 


In terms of this functional , (3-10) implies 


F(p,$ p ,U) 



(3-16) 


and (3-11) implies 

F(p',$ P ,U) > F(p,$ p ,U) for p ? p' . (3-17) 

Finally, since P p is a unique functional of the p^ (r) , (3-15) 
may be written 


F(p,$,U) = 4>i (£,) P^r) + G (p , U) 


(3-18) 


where 


G(p,U) = JdR P P (R){U(R) + ~ In P P (R) 

+ J £ ln(N. !>} 

g i 1 


(3-19) 


is a unique functional of U and the densities p^ (r) . This is 
the function G introduced in equation (3-6) , as can be seen by 
setting $ = $ p in (3-18), and using (3-16). 

As we described after (3-6), since G(p,U) is a unique 
functional of p, we are at liberty to write a formal expansion 
of G(p,U) about its value in the uniform system, specified by 
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the densities p®. We proceed as follows. Consider that for 
0=0, {3-18) reads 

F(p,0,U) = G<p,U) . (3-20) 

Further, since the p? are the densities appropriate to 0 = 0, 
(3-17) yields 

F(p,0,U) > F (p° , 0 , U) for p ^ P° , (3-21) 

Where p° denotes the set of densities p?. From (3-20) and 
(3-21), we see that G(p,U) satisfies 

G ( p , U ) > G(p°,U) for p / p° . (3-22) 

G thus has a stationary point at p?, and an expansion about 

* 

this point must be of the form 
G(p,U) = G (p° ,U) 

+ 5 -Hfd£dr’ fip i (r)6p 1 (r')f ii (r-r') (3-23) 

i j J J J 

+ o (6p) 3 . 


*We note that the relationship between G(p,U) and fj^ is not 
unique. In the canonical ensemble from which we have derived 
these results, the particle numbers Ni are fixed, so that for 
any allowable pj_(r) (which can be obtained by the application 
of an external field) , we must have 

Jdr 6p i (r) = 0 . 

Then on the right side of (3-23) , the addition of any constant 
term to the f^j will not change the value of the integral. 

This difficulty is not important, however, because the physics 
of our problem is not contained in the results of setting k 
rigorously equal to zero (since these results are ensemble 
dependent) , but only in the results in the limit as k 
approaches zero. As the development proceeds, we shall thus 
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It is the functions f^ which play the central role in our 

work. For the case m = 1, the single function f reduces to 

the function f defined in section II. To see this, and make 

use of these functions , we consider now how they relate to the 

response functions x^j • We use (3-23) and the definition of 

_ u . 

the response functions to generate an expansion of F $ m 
powers of the , and compare with the result of a direct 
expansion of (3-1) . 

From (3-18) and (3-16) , we have 

F U p = G (p ,U) + Y2 J" d £ ^(rOp^r) • (3-24) 

If we insert the expansion (3-23) of G(p,U) and make use of 

G(p°,U) = fJ =q , (3-25) 

which follows from (3-20) and (3-16) , we find 

F U p = F$ =0 + £ ♦ i <0) pj + 6p i 

vu i 1 i v k 

(3-26) 

^E5vpp l i««p„(-w*ii„ , y - o ( 6p) 3 . 

Here, we have rewritten the integrals from (3-24) and (3-23) 
as sums in reciprocal space. (4>?(0) denotes the zeroeth 
Fourier coefficient of (£.) • ) We can now pass from an 
expansion in 6p^(r) to an expansion in the <Jk by expanding 
6p i in (3-26) in powers of using (3-4). This yields 

(continued from previous page) ignore further difficulties 
of this nature. This policy will result in several equations 
to which, for precision, one should add the words "plus terms 
independent of r. " 
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f£ = f“ =0+ E+.(0)p? 

+ ¥j + i iS f * m - >Xzi ( - )x I »j t 'y ) 

+ o ( 4>) 1 * 3 . * ( 3 - 27 ) 


Since the functions (r) have now disappeared from the 

formulas, we have dropped the superscripts p. The expansion 
• U 

is now an expansion of F^ in powers of the 

But the expansion in powers of the 4> . may be generated 
without going through the Sp^ expansion and hence without 
introducing the functions f . . . We simply replace the 
potential $ in (3-1) by A4> and expand in powers of A. 


3F 


U 

A$ 


3 A 


JdR *(R)e-^ R > e - eX *< R > 

JdR e -eU(R) e -SAHR) 

£ J d £ (£) (£) > 


(3-28) 


where the single particle densities have been introduced just 
as in passing from (3-12) to (3-13). For the second deriv- 
ative, we formally differentiate the right side of (3-28) and 
make use of the linear response form (3-4) : 




3A‘ 


1 k 

m hT, <fri(-k) + .(k) Xi ^ (k) + o (A) . 

ij v k 1 3 - 


(3-29) 
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TJ 

Then writing a Taylor series for about X - 0, and setting 
X = 1, we have 

F 4 = F 4= 0 + £ 4>i (0) p l 

(3-30) 

+ \ £ h £<M-k)4>. 00x^00 + o (X) 3 . 

ij k J J 

We are now in a position to identify the f^j. Compare 

the result (3-30) to (3-27) , using the facts that these 

expressions must be identical for all potentials 4 , and that 

* 

both x j and f^j are symmetric in k. We must have 

S = -X i3 ty for all i,j. (3-31) 

The structure of (3-31) is made clear if we define the 
matrices F, with elements f^ j , and X, with elements 
Noting from the definition (1-2) of . (k) that ■ - X.^' we 
see that (3-31) reads simply 

X F X = - X (3-32) 

+ 

so that 

F = -X” 1 . (3-33) 


*For Xi j r this symmetry is well known. For the f ^ j , it 
follows because the f^j (r-r 1 ) defined by (3-23) depend only on 
properties of the uniform fluid, hence only on |r-r'|. 

tWe note here another difficulty at k = 0. For (3-33) to hold, 
the matrix X must have an inverse. But the condition 6pi(k=0) 

= 0, which holds for this canonical ensemble, together with 
(3-4), implies that Xi-i (k=0) = 0, so that X” 1 is undefined for 
k = 0. It is, however, defined for all other k, and hence in 
the physically interesting limit as k->0. 
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For a single component fluid, (3-33) reduces to (2-5) , and we 
see that we have recovered the function f whose physics we set 
out to investigate. From (3-33) , the generalization to a 
two-component fluid is immediate: 


-f 


22 


-f 


l ll. 


2 7 X 22 


11 


^ 11 ^ 22 ” ^ 12 ^ 


’ll 22 12' 


x = 12 

^ f f - (f ) ^ 
r ll r 22 lX 12 J 


(3-34) 


Here, we have noted that since x.- ^ / (3-33) implies that 

1 J J 1 

F is also a symmetric matrix. 

We note that the functions f^ bear simple relations to 
the familiar direct correlation functions c. . of Ornstein and 
Zernike (1914) . Comparing (2-5) and (3-34) to the expressions 
in Ashcroft-Lekner and AL-II, we find 


f(k> = -^-{l-p 0 c(k) } 


and 


(3-35) 


’ll 


-d-p^n) 


‘22 


**2 


(l- p 2 c 22> 


‘12 


1 

3 C 12 


This completes the formal development for a general 
system. The key results of this work have been somewhat 
scattered by the development. They are equations (3-6) and 
(3-23), which define G(p,U) and its expansion and thus define 
the problem to be solved for the , and equation (3-33), 
relating f ^ to the The argument to this point essen- 

tially represents a realization for a classical fluid of the 
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case of "almost constant density" studied for the quantum 

electron gas by Hohenberg and Kohn (1964) . 

In our problem, we wish to compare the real liquid to 

some reference liquid, specifically to the hard sphere liquid. 

ref 1 

Dividing the pair potential v into v + v separates U into 
U re f + Then consider G(p,U re ^), which by (3-6), is 

given by 


G(p,O ref ) = F° 


ref 

p,ref 


- E Jdr $?' ref (r)p i (r) . (3-36) 


Here, we have explicitly noted that since must be a 

potential which gives rise to the densities p^ (r) , it will in 

general be a different function in the reference and real 

rs f 

fluids. This function G(p,U ) will have an expansion of the 
form (3-23) with f^ replaced by f^®^, where the f*®^ are in 
turn related through (3-33) to the x^-s • Defining 


G X (P) 


E G(p,U) - G(p,U ref ) 


T ref 


(3-37) 


“ F “p - p ”p,ref * £ Pi(£H^(r)-^' ref (r)} , 


and 


f}. E f . . - ff® f (3-38) 

13 13 13 

it is clear that G^(p) has the expansion 

G 1 (p) = G 1 (p°) + | £ Jdr dr- 6p. (r) 6p. (r* ) f ^ (r-r’ ) 

+ o(6p) 3 


(3-39) 
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Equations (3-37) and (3-39) define the problem to be solved 
for the ff . We note that the relation between f 1 and x has 
already been given for a single component fluid (eq. (2-7)). 
For the two-component fluid, we find from (3-34) 

ref 


ref 

X 11 = ~ ^ 1+ ^ f 22 //f 22 ^ ' x 22 


ref 

*12 ’ ~ ' 




D 


(3-40) 


where 


D = 1 - Y ref f^ - Y «f fl 
X 11 r ll x 22 r 22 


0 ref ,1 
2 x 12 1 12 


k 12 


.ref 


(3-41) 


+ (X^ f xS f - - < f 12> 2 ) • 


Note that in x^/ an Y zeroes of f 2 will be cancelled by the 
(from (3-34)), and that similar reaults hold for the 
other functions. Then the search for singularities of the x's 
is the search for the zeroes of D. 

This completes the development of the basic equations of 
our perturbation theory. The basic equations are (3-37) and 
(3-39) , which define the problem to be solved for the f ^ . , 
and the relation (3-33) with its realizations (2-5) , (3-34) , 

(2-7) and (3-40), (3-41), relating the f's and the x's. 


III-B. Analysis of the Basic Equations 

In this section, we study the role played by the func- 
tions ff j in (3-37) and (3-39) , in order eventually to suggest 



-47 


suitable approximations. Because we shall seek forms for 
G^(p), we shall throughout this section consider that the 
p^(r) are given. 

We must specify the form of U(R). It is to be a sum over 
pair potentials : 


”<*> - | E E V <4 - 4> , (3-42) 

where v^^ is the pair potential acting between a particle of 
type i and a particle of type j , and r^ was defined below 
(3-2). We note the symmetry v.. = v... We note that the 

J J 

restriction £^m is harmless but unnecessary in those terms in 
which 

We shall specialize for a time to the case of a single 
component fluid. 

The structure of our equations takes its simplest form in 
the limit as T-*-0. It is thus instructive to examine this 
limit as an idealization, even though no classical fluid 
exists at T = 0. In this limit, the entropy term in the free 
energy vanishes, leaving only the energy terms <U> and <$>, 
where, as in section I, < > denotes a thermal average. As 
described in appendix A, to evaluate <U> , we must introduce 
the two body distribution function P 2 ) defined above 

(1-7) , while we can evaluate <$> in terms of the single par- 
ticle densities p^r). With these considerations, we may write 


F° - ^-jdr dr' v(r-r ' ) p p (r ,r* ) + jdr 4> p (r)p(r) 


(3-43) 


+ terms which vanish as T-+-0, 
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where P U denotes, just as in the last section, that free 
4> p 

energy determined within a constant by p (r) . The particular 
distribution function p p (r ,r ' ) is thus to be calculated from 
P p (R) and is uniquely determined by p (r) . For the reference 
fluid, we have similarly 
,ref 


5 > 


p ,ref = |f d £ d £’ v ref (r-r*)p^ ref (r,r-) 
+ | dr p (r)<f> p,ref (r) 


(3-44) 


+ terms which vanish as T->0. 


We have explicitly noted that the two body distributions 
determined by p (r) are in general different for reference and 
real fluids. 

From (3-43) and (3-44), we can form G 1 (p) of (3-37). We 

find 

G 1 (p) = |[dr dr 1 v(r-r') p p (r,r’) 

" jjdr dr' v ref (r-r 1 ) p p ' ref (r ,r ’ ) (3-45) 

+ terms which vanish as T->0. 

To accomplish the expansion (3-39), we let ) and 

2T6 f 

P^' (£»£') denote respectively the distribution functions 

determined for the real and reference fluids by the uniform 
density p = p 0 . Then, by (3-45), 

G x (p) = G 1 ( Po ) + |Jdr dr' v (r-r ' ) (p p (r ,r ' ) - p° (r ,r’ ' ) } 

- ^Jdr dr' v ref (r-r ' ) (p p ' ref (r,r' ) - p° ' ref (r ,r ' ) } 
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+ terms which vanish as T+Q. (3-46) 

Comparison with (3-39) then yields 

||dr dr ' v(r-r ' ) (p£ (r,r' ) - p°(r,r')> 

- |Jdr dr 1 v ref (r-r 1 ) {p£ ' ref (r ,r ' ) - p° ,ref (r ,r' ' ) } 

(3-47) 

= |Jdr dr' 6p(r)6p(r , )f 1 (r-r') + o(6p) 3 

+ terms which vanish as T-*-0. 

The meaning of this equation is as follows. If we expand the 
left side in powers of 6p, then in the limit as T-*-0, the 
coefficient of the second order term becomes identical to f^“. 
Then we see that to calculate f 1 (r) at T * 0, we must discover 
how to expand p^ (r ,r ' ) and pP' re f(r,r' ) in powers of 6p. An 
equivalent observation has been of importance in the study of 
the degenerate electron gas, which, at metallic densities, may 
be considered for many purposes to be at zero temperature. 

We turn now to finite temperatures, where we must include 
the entropy term in F. We shall derive two results. The 
first results from the effort to solve the finite T problem by 
the kind of energy consideration which works at T = 0. For 
all T, we define a function f* by setting the left side of 
(3-47) equal to 

jjdr dr' 6p (r) 6p (r ' ) f^(r-r ' ) + o(6p) 3 . (3-48) 

That is, we define f^ to be the coefficient at any T of the 
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second order term in the 6p expansion of the left side of 
* 

(3-47) . Since, as noted above, this coefficient reduces to 

1 

f as T-»-0 , we have the relation 


lim f 1 (r) = lim (r ) . (3-49) 

T+0 T+0 

To discover the relation that exists between f^ and fib at 

E 

finite temperatures, we make use of the thermodynamic identity 


- t(||) = E , 


(3-50) 


where F is the configurational free energy of the system under 
consideration, and E is its potential energy. We explicitly 
note that the usual thermodynamic derivative is taken with the 
external field <j> (r) held fixed. To apply this relation to 
this work, we must instead form thermodynamic derivatives at 
constant density p (r) . A straightforward calculation (sup- 
pressing V and N, which are constant in either case) yields 



(3-51) 


th 

where <f>^ is the k Fourier coefficient of <j> {r ) , and the 
notation <f>‘ indicates that all the are held fixed except 


*To see that this expansion must contain no first order term, 
consider that the general first order term will be of the form 

jdr a (r) 6p (r) . 

But a(r) can depend only on the properties of the uniform 
fluid, so that it must, in reality, be independent of r. Then 
we are left with a term of the form 

ajdr 6p(r) . 

But this is identically zero by the canonical ensemble 
restriction of constant particle number. (See note to (3-23)). 
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the one appearing explicitly in the derivative. But direct 
calculation from (3-1) shows that 




1 

V 


P (“k) . 


(3-52) 


Substituting this result into (3-51) , and rewriting the sum in 
reciprocal space as an integral in real space yields 


(H) p - (I)* + (lf) p fe* P <*>■•<£> * (3 - 53) 

(The superscript p has been added to <j> to make explicit the 
dependence implied in (3-51).) Then, using (3-50), we may 
write 


F - 
We find. 



E 


(lr) p J d £ <f> p (r) p (r) 


(3-54) 


(1 - t(|^) )F U = jjdr dr' v (r-r ' ) p P (r ,r • ) 

' ’ p 

+ (1 - t(|^) ) Jdr <J> P (r) 4> (r) 
P 


(3-55) 


where we have formed the potential energy E just as discussed 

above (3-43). If we now operate with (1-T ) on from 

P _p 

(3-37), and use (3-55), the terms in <j> p and 4> p ' re are seen to 
cancel, and we are left with 


:1 “ T (!t) >V p) = d £' v (£-r’) p p (r,r' ) 


(3-56) 


- |Jdr dr' v ref (r-rMp p ^ ref (r,r') . 
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This implies 

d-T(y ){G 1 (p)-G 1 ( Po )} 

P 

" |Jdr dr' v(r-r' ) (p£ (£/£’) -pj (r,r' ) } (3-57) 

- jjdr dr' v ref (r-r'){p^ ref (r,r*)-p5 ,ref (r,r')>. 

If we now expand both sides of this equation in powers of <5p, 
making use of (3-39) for the left side and (3-48) for the 
right/ and equate the second order terms, we find 


d-T(y )ijdr dr' 6p <r) <Sp (r ' ) f 1 (r-r ' ) 

P 

= |Jdr dr ‘ 6p (r) <5p (r ' ) f* (r-r ' ) . 


(3-58) 


Since this must be true for all <Sp, and since f 1 and fi must 

E 

be symmetric in r, we must have 

(1-T f^fV) = fg (r ) . (3-59) 

(We need no longer form the derivative at constant p, since 
f 1 depends only on the properties of the uniform fluid. ) This 
equation is one of the desired finite T results. It is a 
temperature differential equation to be solved subject to the 
boundary condition that f 1 = at T = 0. It shows the added 
complexity introduced into the problem by the necessary inclu- 
sion of the entropy term in the free energy. 

We can, however, derive a second expression which shows 
that we can regain the conceptual simplicity of the T = 0 
problem if we are content to calculate only to first order 
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in the perturbation v . We make use of the well known tech- 
nique for expanding the free energy in powers of (Zwanzig, 
1954). For an m-component mixture, we have from (3-1) 


N, K.N I j dR e 
1 m J 


-eu(R) -| E (k) $!? (-k> 

* V i V 1 — 1 — 


(3-60) 


where we have written $ (R) in terms of the Fourier components 
of the 4>j_ (£,) and of the density operator defined in (1-3). 
Then, if we also write (3-60) for the reference fluid, we find 

-f^p-F^ref) 


^-SU-^R) -| Z 2$ i (k)*J ,1 (-k^ 


P f ref 


(3-61) 


♦J #1 (r) 


“ <J>? (r)-<fr? /ref (r) , U 1 (R) - U(R)-U ref (R) , (3-62) 


l — T i 


P /ref 


JdR A(R) e -M P ' ref !«) 

fdR e -6U ref (R) e -B* P ' ref (R) 


(3-63) 


In this notation, (3-37) reads 


F $P - F $p,ref = C5 x (p) + p.JkHP'Vk) . (3-64) 


(Note that p^(k) in (3-64) is not an operator.) Now from 
(1-3) , we have 

P ± (k=0) = N i = p i (k=0) . 


(3-65) 
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Inserting (3-64) and (3-65) into (3-61) , we find 


-Bte-^pJ+l ECP i (k) <f>?' X (-k) } 
i k 


,:L (o) /-eu 1 (R) £ £ p. (k> <f>? ' 1 (-k ) 1 

= e> i x 1 /- ^ v i k^O i-i - 


(3-66) 


\* 


1 


ref 


Cancelling the terms in k = 0 yields 


■&{G 1 (p)+i E Ep. (k)*?' 1 ^) 
1 v x k?Q i—i — 




-3U 1 (R) -§ £ £ Pi ( k) <f>?' 1 (" 

^ v i k^O i—i 


k)\ 


(3-67) 


'p , re f 


For U (R) , we have 


ol(R) = *SS v ^4-4> ' 


(3-68) 


where the ^ are the perturbing pair potentials. We intro- 
duce an expansion parameter X by replacing vk by Xvf^ in 
U^(R). We may then formally introduce the X expansions 


G JL (p) = XA(p) + X 2 B(p) + . . . 

4>£ ,;L (k) = Xa i (k) + X 2 b i (k) + ... (k ^ 0) , (3-69) 

U 1 (R) XU 1 (R) . 


If we insert these expansions into (3-67) , and expand all 

around in powers of X, we find (when we note that <j3. (k) > _ 

i — p / rer 

= p^ (k) by definition) 

A(p) = <U 1 (R)> pref . (3-70) 
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This expression can be written in terms of the two body 
distribution functions of the reference fluid. In analogy to 
the notation already introduced for pure fluids , we denote by 
p? j (£ '£ ' ) t ^ e two function, for one particle of type i 

and the other of type j , which is uniquely determined by the 
densities (r) . (We omit the subscript ,, 2 ,, since it is 
cumbersome and, with two indices i,j, redundant.) Then 
(3-70) reads 

A(p) = y£jdr dr' vfj (r-r , )pP' ref (r,r') , (3-71) 

so that, using (3-69), 

G^(p) - G^(p 0 ) 

■ ix^ZJdr dr’ vj . (r-r ' ) { p? ( ref (r ,r ' ) -p?^ ref (r ,r ' ) > 
ij J J J 

+ o(X) 2 . (3-72) 

(The superscripts 11 ref" and " 0 " are introduced here in strict 
analogy to the pure fluid case.) Finally, if we also expand 

i 

fij (r) in powers of X, 

fjj(r) = Aaf.(r) + o(X) 2 , (3-73) 

we can write the expansion of (3-39) as 
G-^(p) — G-j^(Po) 

= | Sjdr dr’ 6p i (r)6p j .(r){XaJ j (r-r l ) + ...} 

+ o(6p) 3 . 


(3-74) 
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Comp a ring the terms of order X between (3-72) and (3-74) shows 


that a^ . 


satisfies 


XI (dr dr’ 

ij J 



v ij ( £"£'HPij ref “ Pij ref (£>£') ) 

d£' 6p ± (r) <5p j <r’ ) (r-r ' ) + o(6p) 3 


(3-75) 


That is, when the left side of (3-75) is expanded in powers of 
6p, the coefficients of the second order terms are the aj^ 
defined by (3-73) , i.e. the first order terms in the expan- 
sions of the f ^ . in powers of the vf . . We note that, as 
x j i D 

promised, the form of this equation is quite like that of 
(3-47) . 

From a calculational standpoint, (3-75) represents 
perhaps the central result. As indicated at the end of sec- 
tion II, the success of the mean field calculation suggests 
that we seek a theory of phase separation by completing the 
low order approximation to f\ Equation (3-75) shows that the 
problem of calculating ff ^ to first order in vf ^ is equivalent 
to calculating P^j re ^ (£#£') t the two body distribution func- 
tion .in the non-uniform reference fluid, to second order in 
the 6p^ (r) . 

From the standpoint of physical approximation, this 
result offers a certain conceptual simplicity. As we noted 
in section II , a more conventional expansion of (k) (hence 
of f^j (k) ) encounters in the first term the three and four 
body distribution functions of the reference fluid. What we 
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have done is to provide an alternate form of the problem. 

Thus, one may either seek to evaluate a^ . by evaluating p?, 

1 3 2 

p|, and p£ in the uniform reference fluid, as required by the 
conventional theory, or by evaluating p^ in the non-uniform 
reference fluid, as required by (3-75). While p° amd p| are 
wholly unknown, we might hope that our extensive knowledge and 
physical understanding of p£ could be used to suggest success- 
ful approximations to p£- This is the approach which we 
ultimately adopt. 

The important results of this section are the "zero tem- 
perature" equation (3-47) , the finite temperature differential 
equation (3-59) (together with the boundary condition (3-49) 
and the definition (3-48) of f^) , and the first order equation 
(3-75). We turn now to the study of approximate solutions to 
these equations. We shall begin by reviewing, in the light of 
these formulas, the matters raised in section II. Then, after 
considering the possibility of a self-consistent approach, we 
turn to the development of the first order formulas on which 
we base our calculations. 



1 1 

IV. Approximate Solutions for the and . 

IV-A. The Earlier Results 

The simplest approximate solutions for f^j are those 
which give, rise to the mean field theories of section II. We 
shall specialize, as in section II, to the case of single 
component systems. 

Consider first the case in which the reference system is 
an ideal gas. For this case, we proposed in section II the 
“simplest mean field theory" specified by (2-12) or 

f 1 (r) = v (r) . (4-1) 

This approximation can be seen from our formalism to be the 
correct first order approximation to f^ (a result which is 
well known from the perturbation series.) Consider the first 
order equation (3-75) . For a pure fluid, it reads 


Jdr dr ' v^r-r') {p£' ref (r ,r 1 > - p° ,ref (r ,r ' ) > 


-/■ 


dr dr ' 6p (r) <$p (r ' ) a^r-r') + o(<$p)^ . 


(4-2) 


When we take the reference system to be an ideal gas, we have 

ref 1 , 

v = 0 , v = v , and 


pP,ref(^^,) _ p (r) p (r ' ) . 


(4-3) 


Thus, the left side of (4-2) becomes 
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where have made use of the canonical ensemble restriction 
(see note to ( 3—2 3 ) ) 

Jdr Sp(r) = o . (4-5) 

Comparing (4-4) to the right side of (4-2) , we conclude 

a 1 (r) = v(r) , (4-6) 

where we have assumed that v(r) = v(-r) (an assumption we 
shall make throughout this work) and required (4-4) and the 
right side of (4-2) to be identical for all 6p (r) . Then 

f X (r) = v(r) + o ( v) 2 , (4-7) 

so that (4-1) is just the first order approximation to f 1 (r) . 

We note that (3-75) allows a ready generalization of the 
simple mean field theory to multi -component systems. For 
v£J £ = 0, (3-75) yields 

a ij<£> “ v ij<£> ' <4 ~ 8) 

(To produce (4-8) from (3-75), we must consider varying 6p i (r) 
independently for each i.) Then to generalize (4-1), we just 
make the first order approximation 

f*j(r) = v^j (r) . (4-9) 

The resulting partial structure factors for the two component 
case are precisely those used by Stroud (1973) in his study of 
phase separation. 

The approximation (4-1) can also be understood in a more 
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sub tie , but instructive, way from our equations before the 
first order approximation is made. Consider first the equa- 
tion (3-47) , which represents the limiting case as T-»-0. For 
the case of a pure fluid and an ideal gas reference system, 
this equation reduces to 


if 


dr dr' v (r-r 1 ) { (r,ir' ) ~ ^ 

= |jdr dr’6p(r)6p(r') f 1 (r-r ’ ) + o(6p) 3 , 


(4-10) 


(where we have ignored the terms which vanish as T-»-0.) Now 
the and p° are in general unknown functions, since they are 
the functions appropriate to the "real" fluid. But because 
the liquid exhibits no long range order, we do know the form 
of these functions when the points designated by r and r' are 
sufficiently far apart. Then the functions just assume the 
form (4-3) : 

p^(r,r' ) = p(r)p(r' ) . (4-11) 

We shall speak loosely of the particles being "correlated" in 
the short range in which p^ (r ,r‘ ) differs appreciably from 
(4-11). Now suppose in (4-10) that the potential v(r) is well 
behaved at small r, and of sufficiently long range that the 
dominant contribution to the integral on the left of (4-10) 
comes from regions beyond the range of correlations, that is, 
from regions where (4-11) is a good approximation. Then we 
might reasonably ignore the short range effects, and replace 
p 2 everywhere by p(r)p(r'). On making this replacement in 
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(4-10) , we recover immediately the mean field form (4-1) . To 
extend this result to finite temperature, we note that the 
left side of (4-10) is just the (implied) left side of (3-48) 
for this case, so that the approximation (4-11) yields immedi- 
ately f3(r) = v(r) . This result is independent of temperature, 
so that solving (3-59) subject to the boundary condition 
(3-49) yields again (4-1). Thus, the mean field theory of 
(4-1) (or (2-12)) emerges at all temperatures from the neglect 
of short range correlations in {£/£.' )• 

Next, consider the case of an arbitrary reference system, 
for which we proposed the mean field form (2-14) , or, more 
precisely , 

f 1 (r) = v 1 (r) . (4-12) 

This approximation results in our formalism when we make the 
first order approximation f^ = a^, and evaluate a^ in (3-75) 
with the approximation 

p£' ref (£,£') = p (r) p (r 1 ) . (4-13) 

In words, for an arbitrary reference system, (4-12) represents 
f^ calculated to first order in v^ with neglect of reference 
system correlations. 

We are now better able to understand the modification 
(2-16) which was introduced to remove the ambiguity of (4-12) 
(or (2-14)) in the special case v re ^ = v HC . Since two parti- 
cles in a hard core fluid can never be found separated by a 
distance less than a (where a is the hard core diameter of 
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D HC 

(1-10)), we must have (r,r_' ) = Or identically, when 

| < a. It is clear that the simplest approximation in 
the spirit of (4-13) which satisfies this condition is 


P 


P / HC 
2 


(r , r 1 ) 


( 0 | r-r ' | < a 

p (r) p (r' ) | r-r '[ > a 


(4-14) 


When a^ is evaluated from (3-75) with this approximation, the 
form (2-16) results immediately from the first order approxi- 
mation f 1 = a 1 . 

We are now better able to understand the significance of 
the correlation effects at liquid densities. At these densi- 
ties, the hard sphere correlations have a range comparable to 
that of v^, so that the dominant contribution to the integral 
on the left of (3-75) (or (4-2)) comes from the region where 
(4-14) is a poor approximation. 

Consider next the limit of close packing, which was 
introduced in section II to illustrate the correlation 
effects. Within the formalism we have developed, the case of 
close packing assumes a special importance, for it represents 
the one case, apart from the case of an ideal gas, for which 
we can write an explicit form for (£/£')• In a close packed 
system, we may induce a non-uniform density p (r) by applying a 
field <j> (r) which has the periodicity of the lattice. But this 
density cannot change the conditional probability that a par- 
ticle will be found at r' given a particle at r. Since this 
probability is given by (r ,r ' ) /p (r) , we have 
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p£(r,r») p®{r,r' ) 

p(r) Po 

or 


(4-15) 


p?(r,r') = — - pf(r,r') , close packing limit, (4-16) 
P O * 

We see immediately that (£*£') has no second order term in 
6p. Though some care is required/ it is no surprise that 
these considerations can be made to yield f*"(k) = 0 when k is 
any reciprocal lattice vector, and hence, the desired result 
X = x • Thus we have achieved one goal we set for this work, 
namely, to handle the close packing limit correctly. 

Return now to the problem at liquid densities, and con- 
sider a more sophisticated approach than is represented by the 
mean field theories discussed so far. Consider the structure 
of equation (4-10) , which holds for the case of an ideal gas 
reference system in the limit as T+0. It states that when the 
expression on the left is expanded in powers of 6p, the func- 
tion f^ will appear as the coefficient of the second order 
term. In order to evaluate the expression on the left, how- 
ever, it is necessary that we already know f^". This follows 
because p°(r,r') (which is the two body function in the uni- 
form "real" fluid) is related through g(r-r' ) to S (k) and 
X(k)> and is thus sufficient to determine f 1 . If a stronger 
statement were true, if the knowledge of f^ were sufficient to 
determine the entire left side of (4-10), we could define a 
self-consistent procedure for determining f^ through (4-10). 
Such a statement is, of course, not true without further 
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approximation, since f"*" says nothing directly about p^Cr/r’ ). 

Similar observations on the problem of the quantum elec- 
tron gas have led Singwi and co-workers (1968, 1970) to sug- 
gest approximate forms for (£,£.')/ for which, in the present 
context, a knowledge of is sufficient. The earliest of 
these was the form 

p£(r,r') = p (r) p (r 1 ) g° (r-r ’ ) , (4-17) 

where g° (r) is defined (as in section I) by g° (r) =p ° (r) / (p 0 ) 2 . 
With such approximations, Singwi et al. were able to define 
self-consistent methods for calculating the electron gas 
response functions. Within the classical formulation repre- 
sented by (4-10) , the particular approximation (4-17) yields 

f X (£) = v(r)g°(r) , (4-18) 

which is to be solved self-consistently for f 1 and g°. We 
note in passing that, by (4-16), the particular approximation 
(4-17) fails in the limit of close packing, so that it is 
highly suspect in application to a dense fluid of hard cores. 

The point we wish to consider here is not the particular 
form of these approximations, but rather the possibility of a 
self-consistent approach to our problem. When we turn to 
finite temperatures, we find that even the simple approxima- 
tion (4-17) yields an exceedingly complicated problem. When 
substituted into (3-48) , it will yield 

4 <£>. = v (r) g° (r) . 


(4-19) 
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This is a temperature dependent result, so that, to find f\ 
we must integrate the equation (3-59) , maintaining self- 
consistency throughout. Because of this complication, we have 
not attempted a self-consistent calculation in our problem, 
despite the success which has accompanied such calculations 
for the zero temperature electron gas. 

We have instead elected to pursue our original plan of a 
hard core plus tail perturbation approach, and focus our at- 
tention on the first order equation (3-75) , in which the ex- 
pression on the left is a function only of v"^ and the proper- 
ties of the reference fluid. 

IV-B. The Solution for a^.. 

We begin with the remark that (3-75) can in fact be 
solved exactly for af ^ in terms of reference system distrib- 
ution functions. This can be accomplished by expanding 
P 2 , ref(r,r l ) in powers of the fields $ p which give rise to 
6p (r) . Though this expansion will contain terms of first 
order in $, we see from the form of (3-75) that the first 
terms surviving the integration on the left will be of second 
order in 4>. It is then a simple matter to use the response 
functions x to produce the second order terms in 6p. Now the 
expansion of p^'^^fr^r') in powers of $> p introduces the 
higher order distribution functions p|' and p°' , so that, 

as might be expected, the result of this approach is simply 
the more conventional form of f^ alluded to in the discussion 
following (3-75) . For the reasons outlined in that discussion, 
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We prefer to work with (3-75) directly, and seek approxima- 
tions for P 2 ,re ^ (r,r' ) • 

This task may be simplified by considering the problem in 
reciprocal space. Choose the densities <$p^(r) to be given by 

6p^(r) = 6p^ cos (k^» r) , i = l,2,...m. (4-20) 

When the right side of (3-75) is evaluated for such densities, 
the full equation takes the form 


Xjjdr dr' vj_. (r-r ' ) { p?? ref (r ,r 1 ) - p? ,j ref (r ,r 1 ) > 

= ?S'Sp i 'Sp>af.{k)(6 k +5 v _k^ + o { «S p ) 3 . 

2 ij 1 ^ 13 ^i'^j ^i' 


(4-21) 


where we have introduced Fourier transforms, and the function 
6. , , has its usual meaning. This result means that, in order 

t K 

to calculate the a^ (k) correctly at some k, we have only to 

evaluate p? ( ref (r,r’ ) for density variations specified, 

3 

through (4-20) , by the same k. (The complexities of (4-21) 

are such that some care must be exercised in order to reach 

this conclusion for a multi-component system. ) Then the 

problem of the af . (k) can be solved independently at each 

^ 3 


wavelength. 

1.) The long wavelength limit. 

Consider first the case of a single component system on 
which we impress a long wavelength density variation. Clear- 
ly, as the wavelength of the density variation becomes infi- 
nite, 6p will go to a constant, and p P ' rex (r,r’ ) can be 
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expanded in a conventional Taylor series in 6p: 
pP' ref (r,r') - p°' ref <£,r') 

(4-22) 

= 6p |— p2' ref (r,r’) + |(<5p) 2 — — p°' ref (r,r') + o(5p) 3 . 
Sp 0 2 2 3Po 2 2 

Now for wavelengths which are finite but still very long com- 
pared to the range of correlation, we can simply evaluate 
P 2 ,re f (r ,r ' ) at the local density. We generalize (4-22) for 
this case by 

pP' ref (r,r’) - P °' ref (r,r’> 

* | {Sp(r)+6p(r’)> |^p°' ref (r,r') (4-23) 

+ i 6p (r) 6p (r 1 ) — — P 2' ref (r,r’) + o(6p) 3 . 

2 3Po 2 2 

Within the range of correlations, by hypothesis, 6p does not 

vary significantly, so that (4-23) is identical to (4-22). 

The particular form has been selected from considerations at 

long range. Beyond the range of correlation, where 

p2' ref (r,r’ ) = p 0 2 , (4-23) reduces to the correct form 

p?' ref (r,r’) = p (r) p (r ’ ) . Substituting from (4-23) into 
^ ' ““ “ 

(3-75) , we read off immediately 

a^ (r ) = \ v L (r) — — p°' ref (r) (4-24) 

apo 2 

3TG f 

Here, we have noted explicitly that, because p^' (r ,r 1 ) 
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the function appropriate to a uniform fluid, it is a function 
only of r-r' . In Fourier transform, (4-24) becomes 

=i|— ^v^k 1 ) — pf ,ref (k-k 1 ) . (4-25) 

^ J { 2 tt ) 3 “ a Po 2 . 

This is the form of a 1 (k) which we expect to be valid for 
small values of k. We shall consider the limits on its region 
of validity more precisely in the next section. 

To extend this long wavelength approximation to the case 
of multi -component systems, we adopt the same basic procedure, 
but the argument is complicated by a lack of symmetry among 
the p?_. in a non-uniform fluid. In particular, we find that, 
in general, pjj (r,r' ) ^ Pjj. (£,£'), that is, that p?_. is not 
symmetric in the indices i,j. For instance, if r and r 1 are 
far enough apart, we have 

Pi j <£»£' (£) Pj (£' ) and Pj i (r,r') = p. (r)p. (r 1 ) , (4-26) 

and these certainly need not be equal. On the other hand, in 
the uniform fluid, the distribution functions are symmetric: 

Pij = Pji (£/£’) • (4-27) 

Now, the long wavelength approximation is to be derived by 
treating the non-uniform fluid as locally uniform. We see 
that we must take care that the approximation does not intro- 
duce the wrong symmetry into the functions for the non-uniform 
fluid. We can do this by defining symmetric distribution 
functions for the non-uniform fluid as well. We define 
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Pij ( £'£ r) = I )+ Pji (£'£') ) • (4-28) 

These functions are explicitly symmetric in i,j, and thus, 
as we shall see, can be treated consistently with a "local 
density" approximation. 

In terms of the p"? j , the left side of (3-75) reads 
Efdrdr' vj\ (r-r‘){F^ ref (£,£') - P?( ref (r ,r ' ) } , (4-29) 

j * X J 

where we have used v^j (r) = Vj^(r) (noted below (3-42)). To 
make the long wavelength approximation, we first expand the 
p"{?. for constant 6p: 

J 


^ ref <r,r-> 


- F?! ref 
13 


(r,r’ ) 


■ ? *01 Ipy "ij" f «£'£'> <4 - 30 > 

+ ? ^ «P/P n p Ij ref (£,£’) ♦ o( 4 p ) 3 


where we have used (4-27) and (4-28) to replace p|j by p?^. 
Just as in the single component case, we proceed to evaluate 
the form of p"?j beyond the correlation range, and then insert 
local densities into (4-29) in a manner consistent with this 
form. (We note that, by (4-27) and (4-28) , both sides of the 
resulting equation will be symmetric in i,j , so that the ap- 
proximation is consistent in this sense.) Substituting the 
resulting approximation into (4-29) and thence into (3-75), 
we find 
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aL(k) = |Zf— vjfk') 

13 2 277 ) 3 im ~ ap i 9 Pj 


o P Jm ref (k-k 1 1 ) , (4-31) 


in complete analogy to (4-25) . 

We pause to note that in the limit as k-*-0 , (4-31) becomes 


a 




v L <y 


8p?3p? 


P 


o /ref 

Jim 


(k) . 


(4-32) 


This limiting expression can in fact be obtained by a short 
argument from a grand ensemble. This argument is given in 
appendix C, and provides a check on our work. 

2.) The problem of short wavelengths. 

In section II , we discussed the difficulty encountered 
by the mean field theory (2-16) , when applied to dense fluids, 
at wavelengths around the first peak of the reference (hard 
sphere) structure factor. There we introduced the ideali- 
zation of a close packed lattice in order to interpret the 
difficulty. As we have seen, the formalism we have developed 
reproduces the expected form in the close packing limit. 
Despite this improvement, however, we still encounter consid- 
erable difficulty, at liquid densities, in formulating a 
successful approximation for values of k near the reference 
peak. Reduced to simplest terms, the difficulty is that, 
while an explicit form (4-16) exists for p 2 ' ref {r,r' ) in the 
close packing limit, we must proceed at liquid densities by 
approximation, and the formulas, for values of k around the 
peak, are exceedingly sensitive to errors in f^. This 
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sensitivity results , of course, from the peak itself, since 
ff. appears in the combination f ^ (and also in the more 
complex quadratic term of (3-41) , for which similar consid- 
erations will hold.) In the single component case, for 
example, X HC (k) for k at the peak is nearly 100 times as large 
as for k = 0. 

The difficulty might be measured against the problem in 
the electron gas. One approach which has been successful in 
the electron gas problem has been to define a short range cor- 
relation function, which, in the spirit of our work at long 
wavelengths, is evaluated at the local density. Singwi et al. 
(1970) take this approach in a paper subsequent to that dis- 
cussed around (4-17) . They define and evaluate at the local 
density a correlation function a which is a 6-function in the 
RPA , and which they expect to be short ranged as well in the 
fully correlated solution. Such an approach seems always to 
fail for the hard sphere liquid, because the hard sphere poten- 
tial, unlike the Coulomb potential, defines a characteristic 
length a which is, unfortunately , of order H&k at liquid 

densities. It has not proven possible to define a correlation 
function whose range is smaller than a, and which thus might 
believably be calculated at the local density when k ~ k peaJ <.* 

We did make an attempt to evaluate at the local density the 
direct correlation function of Ornstein and Zernike (see 
(3-35)), which has a range of order a. (The range is exactly 


a in PY.) Not surprisingly, this calculation gave very poor 
results around the peak, and will not be reported in detail. 
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If an answer is found, we expect that it may arise from 
consideration of the finite size of the particles. To test a 
simple idea, we noted that, at liquid densities, the wave- 
length at is very nearly equal to a particle diameter 

and to the average particle spacing. It thus seems reasonable 
to suppose that at short distances, where the correlations are 
strong, a density variation appropriate to k - k , might 
have effects similar to long wavelength density variations. 

(We are here moving in the direction of the reciprocal lattice 
vector concept appropriate to a crystalline solid. ) Noting 
further that the long wavelength form (4-23) is correct beyond 
the correlation range at all wavelengths, these observations 
suggest that this form (4-23) might have a region of validity 
around the peak as well. As we shall report in the next 
section, the calculations yield interesting, but wholly 
inconclusive, evidence that this might be the case. Perhaps a 
more dilligent pursuit of the short range analogy with a 
crystalline solid may yield a successful approximation. 

We note finally that these difficulties do not have the 
effect of terminating (or invalidating) this effort. It is, 
after all, in the long wavelength region of the spectrum that 
we expect to observe the phase separation instability, and we 
have a clear approximation for this region. Further, as we 
shall see, the correlation effects we have learned to include 
are of considerable importance at long wavelengths as well. 
Thus, while there is still much to be learned, the progress 
which has been made in this effort can be put to good use. 
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As set forth in the opening section/ the program of our 
work is to develop a suitable perturbation approach to the 
hard core plus tail mixture, and to apply this theory to 
liquid metal systems by constructing hard core plus tail 
models of the real potentials- With these last remarks on the 
short wavelength problem, that portion of our effort devoted 
purely to the perturbation development comes to an end. We 
turn now to the application of this work to liquid metal 
systems, to learn about these systems, and, of course, to 
learn more about the perturbation theory itself. 



V. Application to Liquid Metal Systems 
V-A. Systems and Models 

To lay the ground work for this section, we must present 
the form of the effective pair potential and define the 

parameters we shall use to discuss both the metallic systems 
themselves and the hard core plus tail models we shall propose 
for these systems. 

The effective ion-ion pair potentials used in this work 
are those which Ashcroft and Langreth (AL-II) have used to 
discuss liquid metal structure. The potential is derived and 
discussed in AL-II. It consists of the direct Coulomb inter- 
action between a pair of ions, and an indirect interaction 
between one ion and the polarization induced by the other in 
the gas of conduction electrons. In reciprocal space, it 
takes the form 


, T 4ttZ . Z .e 2 

v^(k) = ±-X- 

13 ~ k 2 


{1 + 


v? s 

X 


(k) V P S (k) 


{4iT/k 2 ) 


l £(k) 


- 1 )} 


(5-1) 


Here, we have presented the interaction between an ion of 
species i and an ion of species j. Z^ and Z^ represent the 
respective valences, e is the electronic charge, e (k) is the 
conduction electron dielectric function, and v? s (k) is related 
to the bare (unscreened) pseudopotential, (k) , acting be- 
tween a conduction electron and an ion of type i by 

vj(k) = -Z.e 2 vP S (k) . (5-2) 
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To specify the pseudopotential and dielectric function, we 
follow AL-II exactly. For v? s we choose the empty core poten- 
tial (Ashcroft, 1966) given by 


± r>r i 
r c 


v? s (k) = — cos (kr^) , 

1 k 2 c 


where is a parameter to be specified for each species. We 

have used in this work the values of r^ listed by Ashcroft and 

c J 

Langreth (1967C) . A more complete listing of suggested r 1 

C 

values is given by Cohen and Heine (1970) . (See also Stroud, 


1973) . For the dielectric function, we use 


e (k) = 1 + ■— F(y) , y = , 

y 2 4K F 

F(y) = f (y) {1 - X 2 f (y)/(2y 2 +g) } _1 

f (y) = j + 17 U-y 2 ) ln |x^y | (5 " 4) 

g = (1 + 0.158X 2 ) -1 

X 2 = (ira 0 k F ) _1 , a„ = , k p = (3ir 2 p e ) 1 / 3 , 

me 2 

where p Q is the conduction electron density. These formulas 
completely specify the pair potential when the parameters 
r 1 , Z. , and p are known. The real space form of these poten- 
tials for some systems are calculated and presented in AL-II. 

In calculating the potential for a given system, we shall 
consider that the parameters r^ and are characteristic of 
the species. Beyond the specification of the components, 
then, the parameter from which all else follows is the 
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conduction electron density p . In this work, we replace p 


by the usual dimensionless parameter r , defined by 

5 


S 


4Tra 0 3 ^e 


( 5 - 5 ) 


This parameter is the radius, in atomic units, of a sphere 
with volume equal to the average volume per conduction elec- 
tron. For a pure metal, r may be determined from a know- 

3 

ledge of the density and valence of the element. 

The specification of r for the alloy is more complicated, 

3 

since it will be a function of concentration. The alloy is 
specified by the average densities, p^ and p 2 , of ions of 
species one and two. As indicated in section I, we find it 
convenient to replace these variables by 

Po 


p “ Pi + P 2 


and x = 


P 1 *P 2 


( 5 - 6 ) 


To determine the total ionic density p, we shall make an 
interpolation between the densities p? appropriate to the 
individual components in their pure state. The interpolation, 
used also in AL-II, is 


P (x) 


1-x , x 
J? ri P 


( 5 - 7 ) 


If 3Tg and r* are the values of r g appropriate respectively to 
pure species one and two, the equivalent interpolation in r 


is 


(l-x)Z 1 (r^) 3 + xZ 2 (r^) 3 


s 


(l-x^ + xZ 2 


( 5 - 8 ) 
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This completes the description of the metallic systems and 
their potentials. 

We turn now to the description of the hard core plus tail 
models of these systems. Consider first the hard core refer- 
ence fluids. To specify a single component hard sphere liquid, 
we require, in addition to the particle density p, the hard 
sphere diameter a of the potential (1-10) . We find it conven- 
ient in this description to replace cr with the packing frac- 
tion ri defined (as on page 13) by 


volume occupied bv hard spheres 


total volume 


The convenience of this arrangement lies in the fact that 
classical single-component hard sphere liquids of different 
densities but identical packing fractions are simply scaled 
versions of each other. Specifically, the structure factor of 
any such system, expressed as a function of the dimensionless 
wavenumber variable y = ka, is completely determined by r» 
alone. 

To specify the two-component hard sphere liquid, we 
require, in addition to the densities p-^ and p 2 , two hard 
sphere diameters cr^ and a 2 , in terms of which the reference 
interactions take the form 


r<a i j 5 (ff.j+c.j)/2 


r>0 ij 


For the same scaling reasons applied to the single-component 
case, we replace the diameters and a 2 by n and a, defined 
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by 

volume occupied by hard spheres 

n - 6 ^P 1 o 1 + p 2 o 2 ) - total volume 

(5-11) 

ct = o^/a 2 (0 £ a £ 1, by convention) 

In this case, the partial structure factors, expressed as 
functions of the dimensionless variable y = ko 2 , are uniquely 
determined by x, r) , and a. 

In our calculations for the hard sphere reference sys- 
tems , we have used the structure factors determined from the 
PY approximation. (References were listed in section I.) 
Because the forms of these functions are given by Ashcroft and 
co-workers (1966, AL-I) , and because, for our purposes, we 
need only be aware that they exist, we shall (mercifully) 
neglect to copy them into this document. 

Finally, we consider the description of hard core plus 
tail models for the metallic systems. Once the reference hard 
sphere parameters have been selected, we pass to the hard core 
plus tail model by adding to the potential v HC that part of 
which lies outside the core. That is, in the pure fluid 
case, our model potential will be 


v(r) 


f 

} AL, , 
I. V (r) 


r<cr 

r>c 


(5-12) 


while for the binary case, the model potentials will be 
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v ij <r > 


■ r 

(v“ 


r< °ij 

“< r) r>a ij 


(5-13) 


It is clear that the adequacy of these model potentials will 
depend greatly on the choice of the cr's. This choice is, in 
fact, a central problem in this work, and will be discussed in 
detail later on. 

To use our perturbation theory for the potentials (5-12) 
and (5-13) , we need to specify the perturbation potentials 

. It is clear that we have already made the partial spec- 
1 A L 

ification (r) = v^j (r) for r > . The region r < is 

inaccessible, so that the form of in this region is in 

principle immaterial. Indeed, as we shall indicate, this 
feature is explicitly represented in the formulas we have 
derived. Computational convenience, then, (meaning the use of 
(5-1) through (5-4)) usually dictates v^j (r) = v^(r) in this 
region as well. Then, to study the model potentials (5-12) 
and (5-13) with the theory we have developed, we study the 
perturbation problem defined by 


v raf 


HC 


and 



AL 

v 


(5-14) 


V-B. The Mean Field Calculations 

In this subsection, we present the mean field calcula- 
tions referred to in section II. Consider first the line of 
singularities in the x-T plane. This is to be calculated by 
generalizing the mean field formula (2-16) to the case of a 
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binary system. That generalization may be established in 
precise analogy with the discussion around equation (4-14) , 
and is, not surprisingly, 


f ij ( r ) = {1 - j (r ) } v\. (r) , 


(5-15) 


where 


H. . (r) 
xj 


(l r <o , . 

13 ■ 

( 0 r>a . . 

v n 


(5-16) 


With the choice (5-14), the form of ff .. (k) (i.e. the form in 

reciprocal space) is readily calculated by numerical convo- 
lution of v^(k) with the Fourier transform of H. .. 

To find the singularities in the (k) , we seek the 
zeroes of the denominator D of (3-41) . Expressing the denom- 
inator in terms of the structure factors (using (3-5)), we can 
write 


D - 1 + + xs“fj 2 + 2V5TT5T S«=fJ 2 } 


+ 8 2 P 2 x (1-x) { ) 2 


q HC HC 
b ll b 22 


}{(£^) 2 


(5-17) 


f 1 fl 
E 11 22 


} . 


We studied this denominator for the Li-Na system. The valen- 
ces Z. and Z 2 are both unity, and we use r^ = .561 A and 
— — c 

Na 0 

r c ~ *®^4 A. To simplify the calculation, we have taken both 

the density and the parameters a, n to be independent of tem- 

1 HC 

perature. Then both f^^ and are independent of tempera- 

ture, so that (5-17) is a simple quadratic in 8. We calculate 
by interpolating with (5-8) between the representative 



-81- 



FIGURE 3 

Comparison of experimental phase boundary and MF line of 
singularities in Li-Na. The parameters of the MF calcu- 
lation are specified in the text. 
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values r = 3.333 and r = 4.064. Finally, we choose a and 
s s 

n to be independent of concentration as well. The values used 

in this calculation are a - r" x /r* a - .8201, and n = .456. 

s s 

(This value of r) is that at which Ashcroft and Lekner achieved 
the pure fluid fit near solidification.) 

With these assumptions and the mean field form (5-15) , 
the line of singularities is calculated as follows. The 
calculation of f f ^ (k) is, as noted above, a straightforward 
numerical integration. At each concentration, we calculate 
these functions at k = 0, combine with S??(Q) (from PY) to 
evaluate D, and solve the resulting quadratic for the two 
values of T at which D =* 0. Since lim D = 1, the stability 
condition (D > 0) will be satisfied for values of T greater 
than the larger root. It is thus the larger root that we take 
as the desired singularity temperature. The resulting line of 
singularities in the x-T plane is presented in Fig. 3, and 
compared to the experimental Li-Na phase boundary of Schurmann 
(1971) . It is clear that the curve has the proper shape and 
symmetry, and although the values of T are low, they are 
certainly of the right order of magnitude. 

The second mean field calculation concerns the structure 
factor of pure liquid Na for arbitrary k. In terms of the 
structure factor, (2-7) reads 

S HC (k) 

S(k) = ~p ~ - . (5-18) 

1 + 8pS WC (k) f x (k) 

We evaluate this formula using the f^ ? of (2-16). This ex- 
pression is of the same form as (5-16) and is thus similarly 
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evaluated by numerical convolution in reciprocal space. 

As we have indicated in several places, a good fit to the 
structure of pure Na just above solidification (T = 371°K) can 
be achieved with the (PY) hard sphere structure factor appro- 
priate to n = .456 (Ashcroft and Lekner, 1966? see also Stroud 
and Ashcroft, 1972). This structure factor is displayed in 
Pig. 4a. We note the strong first peak which indicates, as 
discussed in section II, a strong correlated response at 
ku - 7. Then the mean field theory at these wavelengths 
erroneously predicts a large effect from the perturbation v 1 . 
To display this effect, we calculated the expression (5-18) 
in the mean field approximation, taking the reference system 
appropriate to n = .456, and using T = 371°K and v 1 = v^ 
determined by r“ = .884 A and r a = 4.064. In Fig. 4b, we 
have displayed the resulting form of the denominator of 
(5—18). The strong (for these parameters, catastrophic) 
effect at the peak is clear. 

The significance we have attributed to these mean field 
calculations was discussed in section II, so that no further 
discussion is called for here. 

V-C. Pure Fluid Calculations 

In this subsection, we use the theory developed in 
sections III and IV to calculate the structure factors of 
single-component fluids. These calculations are, in 

*This is a bit of an error with respect to our intent to model 

Na. The correct value at T = 371°K is r Na = 4.045. 

s 
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themselves, of less interest than the alloy calculations to 
come# because the success of our theory is limited to the long 
wavelength regions, and it is only for the alloy case that we 
expect any spectacular effects. Nevertheless, these calcula- 
tions do have some intrinsic interest, and they allow us to 
discuss many important features of our theory and calculations 
without the confusion which accompanies the multi -component 
case. 

The plan of these calculations is dictated by the results 
of the last section and is as follows. We must calculate f 
to insert into (5-18) . We adopt the perturbation approach, 
expanding f^ in a series f 3 = a^ + o(v^) 2 , and truncating 
after the first term. We shall calculate a^ using our only 
clear approximation, the form (4-25), and thus limit our 
investigation to small values of k. In the next section, we 
shall adopt this same scheme, suitably generalized, to calcu- 
late the j (k) . With a name derived from the nature of the 
approximation to a^ (actually to ,re ^ (r ,r ' ) ) , we shall refer 
to this structural approximation as the "mean density approx- 
imation" (MDA) . 

To use (4-25) , we must calculate the function 
3 2 HC 

i p“' (k) . Using the relation between p£ and S (k) defined 

by (1-6) and (1-7) , we find 

|^T p5' HC (k) = 2 ( 2 ir ) 3 6 (k ) 

. ,3 c-HC,,. . n 3 2 q HC ... 

+ 2 3Po S + Po 3p^ S ( ~ } * 


(5-19) 
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At this point, we face a rather difficult matter;. Of the two 
parameters, p and a, needed to specify the HC reference sys- 
tem, p is determined by the particular system under investi- 
gation, while a is left to choice. We shall shortly have a 
good deal to say on the subject of the best choice of o. The 
important point for now is that if a is selected by any scheme 
systematically related to the potential v AL , it will depend on 
the density through the r dependence of this potential. A 
question which might be asked is, should the density deriva- 
tives in (5-19) properly include density derivatives of this 
"best" o? This question is actually the first manifestation 
to this point of what is ultimately a fundamental difficulty 
in applying this work to metallic systems, where the notion of 
an ion-ion pair potential depends crucially on . the accounting 
of conduction electron effects. We shall discuss this matter 
in section V-F. For now, however, we proceed from the assump- 
tion, stated in section I, that the ionic structure factors 
are those appropriate to a classical fluid with pair inter- 
actions v^. Then once a potential has been constructed, 

and subsequently modelled by a hard core plus tail, we shall 
investigate the appropriate classical structure without fur- 
ther consideration of the origin of these potentials. Within 
this view, if the parameter a is accepted as a characteristic 
of the potential v AL , it is a given in the problem, and the 
differentiations in (5-19) are performed only with respect to 
the explicit density dependence. 

We can thus use the definition (5-9) of n to replace 
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% ft 

p ■— by Thr~ • On making this change in (5-19) , and substi- 
3Po »T1 


tuting into (4-25) , we find 


a 1 (k) = v X (k) 


(5-20) 

+ i f_^_ s HC (k-k ( ) + a 2 |L- s HC (k-k*)} . 

P °J ( 2 7T ) 3 - 3T1 2 3n 

up 

We note further that since p 2 ( r ) = 0, identically for r < o, 

ft ^ HC 

we must also have — -rr pS ' (r) = 0, identically, for r < a. 

9po 2 2 

Then (remembering that (4-25) is just the Fourier transform of 
(4-24)), we see that a 1 (k) in this approximation is explicitly 
independent of the form of (r) for r < a. Thus, as dis- 
cussed earlier, we make the replacement v = v . Finally, 
making the first order approximation f^" = a^, we have, for the 
mean density approximation to f^. 


f 1 (k) = v^k) 


- fJL v^k'Hnfrr S HC (k-k') + ^ S HC (k-k ’ ) > 

Po J (2tt) 3 - 3T1 2 3n2 


(5-21) 


We note that the n derivatives of the structure factors can be 


given an explicit analytic form in the PY approximation, so 
that the expression in (5-21) is calculated by numerical 
integration of analytic functions. 

We first specialize this calculation to the case of the 
absolute long wavelength limit, k = 0. Our interest in this 
case derives from two sources. First, at k - 0, the expres- 
sion (5-20) for a 1 (0) is exact (for the hard core plus tail 
potential) . Then the only approximations we make in this case 
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are the first order approximation, and, of course, the PY 

HC 

approximation for S (k) . Secondly, at k = 0 , the physical 
meaning of our results can be discussed using the well-known 
compressibility sum rule: 


lim S (k) 
k+0 



(5-22) 


where the expression in braces is the isothermal compressi- 
bility. 

We have calculated S(0) in the MDA for the v AL appropri- 
ate to four alkali metals, Li, Na, K, and Rb, at, in each 
case, the temperature appropriate to melting. To specify each 
metallic system, we need to specify r , r , and Z. The values 
of Z" s listed in Table I are calculated for each element from 

the density at melting (Smithells, 1967). The r value for 

c 

Rb is taken from Al-Il, while the other r values are those 

c 

listed in Ashcroft and Langreth (1967C) . The valence Z is 
unity in each case. For the hard sphere reference system, we 
follow the suggestion of Ashcroft and Lekner (1966) that each 
of these elements at melting should be characterized by the 
same packing fraction. We therefore take the reference fluid 
in each case to be a hard sphere liquid of the appropriate 
density with n = .456. Within the PY approximation, this 
yields S HC {0) = .02396. 


The results are presented in Table I. The first three 

entries after each element are the melting temperature T^, r g 

° 

at melting, and r in A. The next column presents the values 
calculated for f^(0) according to (5-21). We must say a word 
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S HC (0) = .02396 (n = .456) 


Metal 

t m (° K ) 

r s 

r o (A) 

f MDA (0 ^ 

l+S^f 1 

S MDA*°* 

Li 

453 

3.318 

.561 

-.1965 

.6342 

.03778 

Na 

371 

4.045 

. 884 

.0872 

1.133 

.02114 

K 

337 

5.017 

1.132 

-.0731 

. 9200 

.02604 

Rb 

312 

5.371 

1.12 

-.5247 

. 4587 

.05223 




TABLE I 




Results at k 

& 

o 

II 

= T M ' for 

the pure 

alkalis 

in MDA. 
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about the dimensions used for the entries in this column. In 
this, and all subsequent reportings of f^, , and v^*, we 

adopt the convention of reporting not these figures directly, 
but the dimensionless ratio between these figures and the 

n 

quantity 47rZ 2 e 2 / (k Tp ) 2 = ^ ZE p /p 0 . (Here, k Tp is the usual 
Thomas-Ferrai wavenumber, or, in the notation of (5-4), just 
2k^A 2 .) Thus, for example, in the column labelled f^(0), we 
actually quote the values of, say, £(0), where 

f ( °) - e(0). (5-23) 

1 TF 1 

In these terms, the important quantity 8p 0 f^(k) is given by 
2 

Z£(k)T p /T, where T p is the usual Fermi temperature. The 
last two columns present, from these values of f'*', the values 
of the denominator in (5-18) and the resulting S(0). 

We note from the compressibility sum rule (5-22) and the 

HG 

comparison between S(0) and S (0), that the compressibility 
can be either increased or (in the case of Na) decreased by 
the addition of v . Numerically, this reflects the fact that 
the calculated values of f^(0) are of either sign. We note 
that, for these potentials, the mean field form f^ p (0) (which 
is, from (2-16), just /“v^ (r ) r 2 dr ) will not have this proper- 
ty, since the potentials v^ defined here are predominantly 
negative for r > a. (These potentials are displayed in AL-II. 
In relation to the hard sphere diameter defined by r) = .456, 
they each take a form qualitatively like that displayed in 
Fig. 2.) That the correct from of f^(0) should in fact permit 
either sign is apparent from the fact that these negative v^ 
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neve r the less contain both attractive and repulsive regions. 

At sufficiently high density, the nearest neighbor position 
will be placed in the repulsive region of v^, so that a de- 
crease in compressibility is expected on passing from hard 
core to hard core plus tail models. This is, of course, a 
correlation effect, since a nearest neighbor position can only 
be defined for correlated particles. It is instructive to 
note, from AL-II, that the nearest neighbor position in Na is 
indeed just into the repulsive part of v\ while in Rb, it is 
well into the attractive region. 

We turn now to the results at k / 0. In Fig. 5, we 
present the functions f ^ (k) for Li, Na, and K, calculated in 
the MDA (that is, by (5-21)). Note that to present only 
f ^ (k) , we need not present the temperature (though temperature 

will affect the careful selection of r . ) As with the k = 0 

s 

results, we plot these functions in units of 4irZ 2 e 2 / (k TF ) 2 . 

The MDA is a long wavelength (first order) approximation 

for f . We expect it to be valid for wavelengths longer than 

the correlation range in the reference fluid. The correlation 

range should be given roughly by the inverse of the width of 

HC 

the first peak in S (k) , since this peak defines the dominant 
oscillation of the radial distribution function g (r) (see 
equation (1-6)). From consideration of the HC structure 
factor of Fig. 4a, we expect, optimistically, that what we 
might call the range of "a priori validity" of the MDA should 
extend only out to ka ~ 1. The interesting feature of the 
results presented in Fig. 5 is that, in spite of strong 



- 92 - 


Ng : r s = 4.064 
7 ) = 0.456 


Li, and K. 
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variations in the shape and amplitude of the curves , the 
curves all go through zero at a value of ka near the hard 
sphere peak ^ cf p ea ] c = 6.8). This means that for a small 
region around the peak, the MDA gives to S(k) values which are 
in accord with our physical expectation that S (k) should be 
little affected at the peak by v 1 . This region is really very 
narrow. In Na, the structure factor has been driven to 
negative values by ka = 7,0. (Note from (1-5) that the struc- 
ture factor is a positive definite quantity.) Yet the exist- 
ence of such a region in each element seems to suggest that 
this effect may be, as discussed at the end of section IV, a 
manifestation of finite particle size. 

Until now, we have been content to neglect the problem of 
selecting the hard sphere reference parameters in favor of 
illustrating, with simple (but not unreasonable ) values, the 
general features of these theories. Though we shall adopt 
this tactic once more to make a point when we turn to the bi- 
nary alloy case, we must ultimately consider the selection 
these parameters with care. It is convenient to begin the 
discussion of this matter before we leave the simpler pure 
fluid case. 

For the pure fluid, we wish to develop a systematic pro- 
cedure for selecting the "best" a of the model potential 
(5-12), given the potential v^. The aim of this work is to 
be able, with this model, to calculate closely the structure 
appropriate to v^. Consideration of the physics involved 
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suggests that, with a good enough perturbation theory, a could 
be defined through an empirical approach. The work proceeds 
in two stages. We first model the potential v AL with the po- 
tential of (5-12), that is, with a potential which is infinite 
for all r less than some diameter a, and equal to the real po- 
tential for r > a. We then calculate the structure by a per- 
turbation theory in which the perturbation is that part of v AL 
which lies outside the core of diameter a. Consider for a 
moment a sequence of such models for a single potential v AL , 
arranged in order of decreasing a. If we could sum all orders 
of the perturbation series, and hence calculate exactly the 
structure appropriate to each model potential, we should find 
in this sequence that the structure eventually becomes inde- 
pendent of a. This will happen when v^ takes on a value at 
r « c which is more than a few kT above its minimum value. 

Then further decreases in cr will affect the model potential 
only in a region which is virtually never sampled. In reality, 
of course, we can calculate the structure only in an approx- 
imation which is based fundamentally on the smallness of the 
perturbation. With decreasing a, the perturbation becomes 
larger, and the approximation will eventually break down. If, 
however, the approximate theory works reasonably well for a 
perturbation of order one or two times kT, we might observe 
some vestige of the ideal result. That is, there might still 
be a small range of a for which the calculated structure is 
independent of a. This structure would clearly be the struc- 


ture we are after. 
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To search for such a region, we repeated, for many values 
of rj, the S(0) calculation for pure liquid Na described above. 
A region of independence does not appear. In fact, the result 
is quite strongly dependent on n. AtT=371°K and for 
r“ = .884 A and r" = 4.064, the values of S (0 ) in the MDA 

G S 

vary almost linearly between .0335 at n = .400 and .0140 at 
0 = .500. The result of this calculation is plotted in 
Fig. 6. This result provides a revealing and not too encour- 
aging look at the adequacy of the approach. 

In any case, this result means that the best we can do is 
try to devise an external criterion for selecting a which will 
bring the calculated model structure closest to the structure 
appropriate to It is clear that the essence of that 

criterion will be to select the largest hard sphere diameter 
consistent with a reasonable modelling of the potential. This 
consideration suggests an approach in the spirit of that taken 

by Ashcroft and Langreth (AL-II) , who selected a as that point 

AL 

in the core region of v (r) at which the potential has risen 
^ kT above its minimum. Such an approach seems to contain the 
essential physics of our problem. In more recent work (e.g. 
Stroud and Ashcroft, 1972) a diameter c was selected for a 
different problem by a variational approach. This approach is 
based on the result (Lukes and Jones, 1968) 

F < F HC + <U> HC , (5-24) 

where F is the free energy of the fluid with internal poten- 
HC 

tial U(R), F is the free energy of some hard core fluid, and 
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<U:> is the thermal average of U(R) taken with the hard core 
HC 

distribution functions. For the case of a pure fluid with 

AL 

pair interaction v , (5-24) reads 

F < F HC + Jjdr dr 1 v AL (r-r , )p»' HC (r,r' ) . (5-25) 

In an effort to calculate the liquid free energy, the right 
side of this inequality can be minimized with respect to the 
hard core diameter o. Though such an approach may seem to 

bear little direct relation to our problem, it, in fact, 

3 

appears to contain essentially the same physics as the j kT 

rule, for it gives very nearly the same results for both the 

magnitude and temperature dependence of a. The advantage of 
the variational approach is that the right side of (5-25) can 
be simply expressed as an integral in reciprocal space and the 
whole procedure automated to a great extent. We make use of 
the generalization of this procedure for binary alloys in the 
next section. 

This completes the discussion of our results for pure 
fluids. 

V-D. Binary Alloys I: The k = 0 Results and Phase Separation. 

In this subsection, we apply the theory of sections III 
and IV to calculate for binary alloys the line of singulari- 
ties in the x-T plane. The elements of these calculations are 
the same as those of the mean field calculation presented in 
V-B. We calculate and search for the highest temperature 

at which D of (5-17) equals zero. In the present work, 
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however, we shall replace the mean field approximation to f^\ 
by the mean density approximation, suitably generalized to the 
alloy case. As in the case of the pure fluid calculations, 
the MDA is exact in first order at k = 0, so that, once the 
hard core plus tail potentials are specified, the only approx- 
imations in the present section are the first order approxima- 
tion and, as usual, the use of the PY approximation for the 


S^Oc) . 
13 v -' 


To construct the MDA for the alloy case, we consider the 
long wavelength approximation (4-31) to af ^ (k) , and proceed as 

we did for the pure fluid case. We can express the deriva- 

£ ^ HC trp 

tives g -p og ' p ' o pj^ (k) in terms of the ;T (k) , the reference 

structure factors. We differentiate only the explicit density 

1 71 T 

dependence, make the replacement v. . = vvr because the form of 

2. 3 i] 

(4-31) is observed to be independent of v^ . (r) for r < ck . , 


-1 _ _1 


and finally make the first order approximation f^ ^ = a^ ^ . The 
result is a set of expressions for f^ (k) which are analogous 
to (5-21) but more complicated. These expressions are pre- 
sented in appendix D. They form the MDA to f^ j . We note that 
for this case as well, the density derivatives may be given an 
explicit form in the PY approximation, so that the calculation 
of af . , though by now a lengthy process, is in principle just 

■** J 

a numerical integration over analytic functions. 

For the first calculation, we sought a comparison between 
the mean field and mean density approximations. To do this. 


we use the MDA in a calculation which is otherwise identical 
to the mean field calculation of section V-B. That is, we set 
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r - 4.064, r = 3.333, independent of temperature, and 
a = .8201, ri = .456, independent of temperature and concentra- 
tion, solve the resulting quadratic, D = 0, for two roots T, 
and plot the higher root. The results are plotted in Fig. 7, 
where, for comparison, we have repeated from Fig. 3 the exper- 
imental Li-Na phase boundary and the mean field results. We 
can see that the asymmetry has been maintained, and that the 
curve has been shifted upward by about 100° in passing from 
mean field to MDA. That the curve has been shifted towards 
the experimental phase boundary is probably not important, 
since the selection of a and ri has been quite arbitrary. This 
calculation shows once again that, although we originally con- 
sidered correlation effects at the peak, they are of impor- 


tance at long wavelengths as well. Also, the degree of suc- 
cess achieved by this calculation tends to support our sugges- 
tion that the lesson of the mean field success is that we 


should seek a consistent low order approximation to f . 


To complete this comparison, we present in Table II the 


numerical results of these two calculations. Within our con- 


vention for reporting potentials in k-space, we list, for each 
concentration, the values of v^{0), ff . (0) from mean field, 
and f^j (0) from MDA. We also list the singularity tempera- 
tures for the two calculations. We observe that in passing 
from mean field to MDA, the f^ (0) are substantially changed, 
while their arrangement, in which f^ 2 lies roughly midway be- 
tween and f 22 t is essentially unchanged. Perhaps the 

second of these observations explains why these substantial 
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X 

M 

v^ L 22 
\12 

1 R 

f «4 2 ) 

f MI) 

V<° K) 

t mda ( ° k > 

000 

1.223 
2.455 
1. 841 

-.138 

-.127 

-.133 


- 


100 

1.190 

2.385 

1.787 

-.150 

-.136 

-.144 

-.202 

-.059 

-.168 

. 265 

349 

200 

1.156 
2.314 
1. 736 

-.161 

-.145 

-.155 

-.204 

-.028 

-.160 

351 

434 

300 

1.124 

2.256 

1.691 

-.172 

-.153 

-.165 

-.201 

.005 

-.152 

394 

490 

400 

1.093 
2.201 
1. 648 

-.182 

-.162 

-.175 

-.213 

.033 

-.145 

406 

514 

500 

1.065 

2.151 

1.609 

-.193 

-.170 

-.185 

-.218 
. 055 
-.141 

396 

511 

600 

1.037 
2.104 
1. 571 

-.203 

-.178 

-.195 

-.223 

.071 

-.138 

369 

486 

700 

1.012 

2.059 

1.536 

-.214 

-.186 

-.204 

-.230 

.081 

-.138 

327 

436 

800 

.987 

2.017 

1.503 

-.224 

-.194 

-.214 

-.236 

.085 

-.140 

272 

358 

900 

.964 

1.978 

1.472 

-.234 

-.202 

-.223 

-.244 

.085 

-.143 

205 

243 

000 

.935 
1.941 
1. 438 

-.244 

-.210 

-.232 





TABLE II 

Comparison of MF and MDA calculations of f^j(O) Li“Na. 
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changes in the ff . give rise to only moderate changes in the 

J 

singularity temperatures. 

We must now turn to the problem of the careful selection 
of the reference fluid parameters a and n- The elements of 
the problem here are the same as in the pure fluid case, so we 
to straight away to the variational technique. For the binary 
system, the realization of (5-24) is 

F < F HC + i £fdr dr’ (r-r ' ) p? ' HC (r ,r ’ ) . (5-26) 

ijj J J 

The integral on the right side of (5-26) can be written as an 

integral in reciprocal space over v^j 1 (k) and the partial 

HC • HP 

structure factors (k) . The hard core free energy F is 

known for mixtures in the PY approximation (Lebowitz and 
Rowlinson , 1964). The plan is to select the best values of 
a and n by minimizing the right side of (5-26) at each concen- 
tration and temperature. 

We have applied this variational technique to the calcu- 
lation of the line of singul ari ties of the Li-Na s ystem. The 
plan of the calculation is as follows. Given the concentra- 
tion and a temperature T, we calculate the density of the 
alloy by interpolating with (5-7) between the published den- 
sities of Li and Na at T (Smithells, 1967). In this way, we 
make some account of the thermal expansion of the system. 

From the density, the three potentials v^ 1 (k) are determined 

x j] 

as specified in V-A, and, in turn, the "best" hard sphere pa- 
rameters, a and ti, are determined by the variational tech- 
nique. With the potentials, temperature, and a and the 
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denominator D of (5-17) can be calculated in the MDA. The 

calculation is then repeated for several values of T, to seek 

* 

out that T for which D = 0. This defines an unambiguous 
procedure for calculating the line of singularities. 

We have performed this calculation for three concentra- 
tions, .3, .4, and .5, of Na in Li, using (as always) 

Na 0 Li ° 

r = .884 A and r = .561 A. The zeroes of D are found to 
c c 

lie quite high, at T = lOOflV, 1040°K, and 975°K respectively. 
(The experimental phase separation temperatures are 574°K, 
576°K, and 560°K respectively (Schurmann, 1971).) The value 
calculated for a by the variational technique seems to be 
almost independent of concentration and temperature at .77. 

The values of r) appropriate to the three points (x,T) pre- 
sented here are .274, .271, and .284 respectively. These are 
not unreasonable values, since, in relation to the solidifica- 
tion value - .45 at 400°K for Na, they represent an average 
change of -.03 for each 100°K, precisely as suggested by Al-II 
(though carried to temperatures well outside the range for 
which that suggestion was intended. ) A more complete speci- 
fication of the densities and the calculated values of a and 
n is presented in Table III. 

To round out the picture of the a and t\ dependence of our 
results, we performed a calculation of a different sort, to 
discover what values of a and r\ are required to get the 
"right" temperatures. Of course, we can't fit two parameters 

*For temperature dependent density and hard sphere parameters, 
the denominator is no longer a simple quadratic in T. 
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X 

T(°K) 

Li 

r 

s 

Na 

r 

s 

a 

n 

D(T) 

. 300 

950 

3.42 

4.28 

.770 

.280 

-.033 


1000 

3.42 

4.30 

.768 

.274 

.005 


1050 

3.44 

4.32 

.772 

.266 

.029 

. 400 

573 

3.34 

4.12 

.771 

.351 

-.608 


773 

3.38 

4.20 

. 762 

.315 

-.233 


900 

3.40 

4.25 

.770 

.292 

-.095 


1000 

3. 42 

4.30 

. 773 

.277 

-.020 


1040 

4.43 

4.32 

.771 

.271 

.005 

.500 

950 

3.42 

4.28 

.774 

.288 

-.015 


975 

3.42 

4.29 

.772 

.284 

.000 


1050 

3. 44 

4.32 

.775 

.274 

.048 


TABLE III 

Results of the variational calculation for Li-Na. Values for 
r^ 1 and r^ a are calculated from published density data 
(Smithells, 1967) , a and n are determined by the variational 
technique, and D(T) (at k = 0) is calculated in the MDA. The 
points presented were calculated in the search for T at which 
D (T) = 0. 
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to a single data point, but, for each concentration, there 
presumably exists a locus of points in the a-n plane for which 
D = 0 at the desired temperature. Then the scheme of this 
calculation is as follows. Given the concentration x (we 
choose x - .4, near the critical value, so that the line of 
singularities lies close to the phase boundary) , the observed 
phase separation temperature is defined (575°K from Schumann, 
1971). Then the density and three potentials are calculated 
for this x and T as described above , and the value of D inves- 
tigated as a function of a and n . In Fig. 8, we present the 
a-ri plane for T = 575°K and x = .4. The boundary in this 
plane divides the region where D is greater than zero (above 
the boundary) from the region where D is less than zero. The 
boundary itself is then the desired locus of points at which 
D = 0. For comparison, the variational values assigned to a 
and ri at this x and T are a = .772 and n = .35. 

The range of what we might consider reasonable values for 
a and n is determined primarily by our physical feeling for 
the packing fraction r|* We do not want to select ri much above 
" .47, since the experience in pure fluids suggests that the 
system would tend to crystallize. Since 575°K is about 125° 
above the melting point of Li (at x = .4, the mixture is Li 
rich), the experience of AL-II might suggest an ri of .41 to 
.42. (The lower value obtained from the variational calcula- 
tion presumably reflects both the presence of Na and the 
effects of thermal expansion.) For these reasons, we have 
chosen to exhibit in Fig. 8 the region between n = *35 and 
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0.3S 0.40 0.44 0.48 


1 

FIGURE 8 

The a-ri plane at T = 575 °K, x = .4/ in Li-Na. 
The denominator D is greater than zero above 
the line, and less than zero below. 



-107- 


n = .46. 

From these results, the selection of a and n emerges as 
one of the central difficulties in our work. It is clear that 
the temperature at which D = 0 is quite sensitive to the val- 
ues of these parameters. The variational, calculation, which 
probably has the best a priori justification, has yielded tem- 
peratures which are far too high. We have seen that the 
"correct" temperature can be obtained with "reasonable" values 
for a and t) but the physical meaning of those values, indeed 
the physical meaning of such a fitting procedure, is unclear. 
The confusion at this point suggests that, while we may ulti- 
mately have learned something of the form of the structure 
factors near phase separation, considerably more work on the 
modelling of the systems {and on the theory itself) may be 
needed before we can predict this transition in real systems. 
Our difficulties here may also derive from the special problem 
of the metallic pair potential, to which we shall turn briefly 
in section V-F. 

V-E. Binary Alloys II j The Results for Finite k. 

We turn now to our results in the MDA for binary alloys 
at finite k. In the previous section, we identified points at 
which 


lim D(k) =0 , (5-27) 

k+0 

where D(k) is the denominator of {5-17). In this section, we 
examine one such point in detail, to discover the manner in 
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which D (k) goes to zero in the MDA as k-*0 , and to observe the 
MDA structure factors themselves. Ultimately, we intend such 
work to form the basis of a resistivity calculation, and we 
conclude this section by suggesting a program for such a 
calculation. 

We have chosen the point defined by x = . 4 , T = 575°K, 
a = .825, and n = .430. The values of a and n were estimated 

from earlier results to give D(0) = 0 in Li-Na at this concen- 

tration and temperature. These estimates are in fact slightly 
in error, so that the value of D ( 0 ) is -.0013. In Table IV, 
we list values of y = k/2k p and the corresponding values of 
D(y) in the MDA. These values of y were selected in another 

unit system, and turn out in the present units to be multiples 

of .01692, a fact which we have indicated in the second column 
of Table IV. In these units, the first peak of the structure 
factor occurs at about y = 1.2, so that the region in which 
the MDA is most to be trusted extends out to roughly y = .2. 
When these values of D are plotted vs. y on a log-log plot, 
they are observed to follow quite well a straight line of 
slope two for the small values of y, but move away from this 
line for the larger values (see Fig. 9) . This suggests that 
the singularity may be of the sort 

D « a + bk 2 . (5-28) 

This type of singularity appears in more conventional mean 
field theories. Because of the error in a and n, we would ex- 
pect for these data to set a = -.0013. This suggests that we 
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y (=k/2k F ) 

y/. 01692 

D(y) 

D(y)-D(l 

0.00 

0 

-.0013 

.0000 

.03384 

2 

.0102 

.0115 

.05076 

3 

.0245 

.0258 

.06768 

4 

. 0442 

.0455 

.08460 

5 

. 0690 

.0703 

.10152 

6 

. 0985 

.0998 

.11844 

7 

.1322 

.1355 

.13536 

8 

. 1678 

.1691 

.16920 

10 

.2526 

.2539 

.20304 

12 

.3418 

.3431 

.25380 

15 

.4748 

. 4761 

.33840 

20 

.6579 

. 6592 

.50760 

30 

.7951 

.7964 

.67680 

40 

.8709 

. 8722 



TABLE 

IV 


Results in MDA for D(y) at y ^ 0. Calculated for Li-Na at 
= .4, T = 575°K , a = .825, and rj = . 


x 


430 



D(y) 
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Log-log plot of D(y) in MDA for Li-Na at x = .4, T = 575°K, 
a = .825 and n = .430. The straight line has a slope of two. 



D(y) - D ( 0) 


2.0 



y (= k/2k p ) 

FIGURE 10 

Log- log plot of D(y) - D(0) in MDA for Li-Ma at x = .4, T 
T = 575°K, a = .825 and n * .430. The straight line has a 


slope of two. 
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construct a log-log plot of (D + .0013) vs. y. This has been 
done, and is presented in Fig. 10. The straight line in this 
figure also has a slope of two, and now appears more properly 
as an asymptote. This observation lends support to the pro- 
posed form (5-28) , viewed as a limit. 

There are two points to be made about these results. 

First, they coincide with the physical expectation that the 

singularity at k = 0 is the first (highest temperature) of the 

possible long wavelength singularities. This must be true on 

physical grounds, since any singularity will render the system 

unstable, while real alloy systems of critical concentration 

are stable down to the critical temperature, at which point 

the structure factors diverge at k - 0. Secondly, while we 

. -2 

can present no direct evidence for a k singularity, we must 

-4 

have a singularity which occurs less rapidly than k in order 
that the resistivity integral (1-1) converge at the lower 
limit. That our long wavelength approximation conforms to 
both these physical expectations, despite the complexity of 
the binary alloy formulas, is a gratifying result. 

In Figure 11, we present the structure factors (k) 
calculated at this point in the MDA. They have been formed 
using (3-40) for the x^j and the relation (3-5) between x^j 
and In the same figure, we have presented for comparison 

the (PY) structure factors of the underlying hard core refer- 
ence mixture (the mixture defined by x = .4, a = .825, and 
n = .430). The MDA structure factors appear to join rather 
smoothly with the HC structure factors at point roughly midway 
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between the origin and the first peaks. This fact is in ac- 
cord with our physical expectation for these systems, and, as 
discussed below, might be put to good use in a resistivity 
calculation. However, the MDA is expected to be valid only 
for small k (y < .2) so that the ultimate meaning of this ob- 
servation is unclear. 

This completes the report on the calculations we have 
performed for metallic systems. Although we have stopped 
short of a resistance calculation, we can at least indicate 
the shape that such a calculation may take. We are interested 
in shedding some light on the conjecture, discussed in section 
I, that the resistivity anomaly observed by Schurmann and 
Parks in Li-Na is a long wavelength effect. As indicated in 
(1-1) , what we need for a resistance calculation is the form 
of the (k) for 0 _< k £ 2k F * What we have is an approxima- 
tion valid for small k (0 _< k < 2k F /7) and the physical idea 
that the structure factors (k) should remain substantially 
unchanged around the peaks as we pass from hard core to hard 
core plus tail model potentials. It is immediately clear that 
we are in no position to do any serious quantitative calcula- 
tion. Yet we might still imagine a calculation which could 
lend qualitative support to the long wavelength hypothesis. 

We might, for instance, simply calculate the contribution to 
the resistivity from that part of the integral in (1-1) for 
which the MDA is valid. More completely, we might calculate 
the structure in the MDA at small wavelengths, and join the 
results smoothly to the HC structure factors around the peaks. 
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thus incorporating what we expect will be the essential fea- 
tures of (k) near phase separation. (As noted above, the 
MDA, somewhat inexplicably, seems to be doing this job for 


us. ) 


The trick in any such calculation, of course, is going to 


be to get the T-dependence right (since we are looking for a 
temperature dependent effect. ) This involves understanding 
both the "normal" (that is, linear) and the anomalous parts of 
the resistivity. We might, however, simplify a first calcula- 
tion by supposing that the "normal" part is ultimately attrib- 
utable to the hard cores, that is, could be dealt with by a HC 
structural model of sufficiently subtle temperature dependence. 
(In view of the successful resistivity calculations of AL-II, 
this may not be an unreasonable assumption.) Under such an 
assumption, the plan of a first simple calculation would be to 
ignore the temperature dependence of the density and the 
"best" hard sphere model, and study only the remaining T de- 
pendence of the difference P^g C “ Pei ec * We Pl an to under- 
take such a calculation in the near future. 


V-F. The Metallic Pair Potential 

In these calculations for metallic liquids, we have sev- 
eral times encountered the conceptual difficulty arising from 
the density dependence of the pair potential. We have also 
encountered a great deal of difficulty in getting our calcu- 
lated results to agree in more than general form with the ex- 
perimental results for the real systems. It seems to us in 
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retrospect that these difficulties may be traceable to an in- 
correct assumption made at the outset, namely, that the ionic 

structure is precisely that appropriate to (classical) parti- 

AL 

cles interacting via the pair potential v . We have not un- 
dertaken a detailed and careful investigation of this problem. 
Yet it seems to us that there is , within the tradition which 
gave rise to the effective potential v AL , a certain contradic- 
tion in usage which suggests that the potential v AL is incor- 
rectly used to calculate the structure factors at small k. 
Whether this is in fact the case is a difficult (and still 
open) question. 

The difficulty appears in the following way. Our calcu- 
lations are based on an approximate accounting of energy in 
simple metallic systems (Ashcroft and Langreth, 1967C, and 
Harrison, 1966). For this discussion, the essential features 
of this approximation are contained in the form 

E (p 0 ,R) = E 1 (p 0 ) + E 2 (p 0 ,R) . (5-29) 

The terms belonging to E^(p 0 ) depend only on the average dens- 
ity p 0 = N/V. The terms belonging to E 2 (p 0 ,R) depend on the 
explicit ionic positions R as well as on p 0 . When is cal- 

culated within a linear screening approximation, the result 
can be cast as a sum over an effective ion-ion pair potential. 

With the choices of dielectric function and pseudopotential 

AL 

specified in V-A, this potential becomes v 

Without going further into the origins and form of this 
expression, we can make the important points here by consid- 
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ering and contrasting how it has been used in the calculation 
of the compressibility and structure of simple metals. In 
each of these calculations , we are considering the effects of 
density variations. To calculate the compressibility, we con- 
sider variations in the average density p 0 . A calculation of 
S(k), however it proceeds, must be equivalent to considering 
spatially varying densities, since S(k) is the density linear 
response function. Viewed in this way, the two calculations 
differ only in that the appropriate density variations are of 
infinite wavelength in the first, and of finite wavelength in 
the second. Yet the customary use of (5-29) in these calcu- 
lations requires a fundamental distinction between these two 
cases, hence, between infinite and finite wavelength density 
variations. Good values of the compressibility are obtained 
if one considers all of the p 0 dependence in (5-29). On the 
other hand, the structure factor is calculated ignoring E^(p 0 ) 
entirely, and also ignoring the density dependence of E 2 (p 0 /R)* 
that is, of v A k. (Of course, p c = N/V is constant in the ca- 
nonical ensemble. Our question will be, is this always the 
correct density to use in calculating the energy?) This re- 
duces the structural problem to the problem of a classical 
liquid with pair interaction v AL , and, as discussed in I-D, 
places the hard sphere success on a good foundation. 

The point we wish to make is that, whatever its successes, 
this procedure rests on too sharp a distinction between finite 
and infinite wavelength. Indeed, for long wavelengths, the 
distinction is blurred. This must be so because the correla- 
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tions (both ionic and electronic) in the liquid are of only 
finite range. Then a density variation with a wavelength much 
longer than the range of correlations must have effects which 
are locally indistinguishable from those of a change in p 0 
which gives the same local density. With this blurring of the 
distinction between finite and infinite wavelengths comes an 
ambiguity in the structure calculation at long wavelength. In 
our calculations, we have followed the traditional method of 
ignoring the p„ dependence and treating v^ as a valid pair 
potential at all finite wavelengths. Yet it seems from this 
discussion as if, at long wavelengths, we must discover how to 
include the effects of E^(p„) and the p 0 dependence of v AL , 
presumably by evaluating these somehow at the local density. 

It is not clear which view (if either) is correct. 

We remark that this discussion is confined to the consid- 
eration of long wavelengths. With respect to these matters, 
the existing short wavelength calculations seem to be less in 
doubt. A density variation of wavelength much shorter that 
the range of ionic or electronic correlations will look over 
this range nothing at all like a change in p D . We expect the 
existing pair potential approximation to be valid for such 
wavelengths. These wavelengths will include the first peak in 
S(k). Then, in particular, the success of the hard sphere 
model, which means primarily the ability of this model to re- 
produce the experimental first peak, is probably still attrib- 
utable to the form of v AL . 

Clearly, a careful and systematic investigation of this 
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matter is required before we can remove the ambiguity at long 
wavelengths, and define a consistent procedure. We mention 
that efforts are being made (Watabe and Hasegawa, 1972A,B) to 
treat the structural and thermodynamic properties of pure 
metallic liquids by considering these systems from the outset 
as binary systems composed of ions and electrons. We suggest 
here that one fruitful line of study for the particular prob- 
lem at long wavelength (and in the present tradition) might be 
to consider calculating the energy of the metal in a perturba- 
tion theory that begins, not with the uniform electron gas, 
but with the non-uniform electron gas studied by Hohenberg and 
Kohn (1964) and Kohn and Sham (1965). 

We note finally that if careful investigation of the long 
wavelength problem proves our calculations to be wrong, they 
may still be viewed as model calculations, which illustrate 
the form of the structure factors near phase separation (and 
thus may still give us a qualitative understanding of the 
resistive anomaly). The developments of sections III and IV, 
although carried out with metallic systems in mind, stand 
quite independently of these systems within the tradition of 
classical liquid theory. 



VI. Discussion/ Conclusions/ and Suggestions for Further 
Work 

This work has been divided into two parts: the develop- 
ment of the perturbation theory of the hard core plus tail 
mixture, and the application to the metallic liquids and 
alloys. Although some items to be discussed relate to results 
from both parts, we shall adopt this same division in this 
last section. 

VI-A. The Perturbation Theory 

We have cast our study of structure in terms of the func- 
tions f^j (k) , which we may take as defined by (3-33). These 
functions are closely related to the more familiar direct cor- 
relation functions of Ornstein and Zernike (1914) (see equa- 
tion (3-35)). Perhaps our most important result for these 
functions is the equation (3-75) , which shows that to calcu- 
late these functions to first order in the perturbing poten- 
1 o ref 

tials v. we must calculate p . ' (*>£') r the two body func- 

tion in the non-uniform reference fluid, to second order in 
6p^(r) . This observation leads directly to the long wave- 
length approximation (4-31). As indicated at the close of 
section IV, a successful approximation for short wavelengths 
may possibly be developed by a carefully chosen analogy to 
crystalline solids . 

In this effort, we have produced a theory which retains 
many of the virtues of the ''simplest" mean field theory of the 
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phase transition itself (that is, the theory defined by (4-9) 
and used by Stroud (1973) in his study of phase separation) , 
but which is, in principle at least, a reasonable theory of 
structure (which the simplest mean field theory is not.) It 
incorporates, in principle, the hard sphere successes, and is 
in the spirit of the most successful of the present theories 
of simple liquids, the perturbation theories on a HC reference 
system. 

We note here that, while our effort has concentrated on 
the structure factors, we have actually been doing thermo- 
dynamic perturbation theory as well. To make this point, we 
consider the form of the Zwanzig expansion of the free energy 
P (in the absence of external fields.) Writing 
U(R) = U ref + U^(R), where U 1 (R) is defined by (3-68), spe- 
cializing to pure fluids, and introducing an expansion param- 
eter X by replacing v 1 by Xv 1 , we find quickly (from (3-1) 
with $(R) =0) 


3F- 


3X 


= <U 1 (R)> X = ^jdr v 1 (r)p“ /X (r) , 


(6-1) 


where < >. and p°' X (r) are evaluated from the distribution 
A 2 ~~ 

function determined by U(R) = U re ^(R) + XU^(R). Then 


3 * F x 

ax 2 


Ik v 1 




(6-2) 


Using (6-1) and (6-2) in a X expansion of about X = 0 and 
setting X = 1, and writing p^'^ -0 ^) ( = ^ xe; ^ (r) ) in terms 
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of S ref (k) by (1-6) and (1-7), we find eventually 


F = F 




C dk 

ref + N|a. { v 1 (o) + i v 1 (k) (S ref (k)-1)} 

2 po ' (2 tt) 3 “ 


1 fdk ! 

£ v ^ 

Po ^(2tt) 3 


(k) (S ref (k)) 2 (-6 Po a 1 (k)) > (6-3) 


+ o ( v 1 ) 2 . 


where a^ (as defined in (3-73)) is the first order term in the 
\ expansion of f 1 . That is, the second order term in the 

3T0 f 1 

Zwanzig expansion of (2-2) may be calculated knowing S , v , 
and a^. At the present level of our work, of course, in which 
we have only a long wavelength approximation to a^ (k) , we are 
not yet able to perform a sensible approximate calculation of 
this second order term in F. Perhaps with better approxima- 

D 270 f 

tions to this situation will change. We note that this 

second order term is, of course, the term studied by Barker 
and Henderson (1967A, 1967B) . In this sense, our MDA appears 
similar in spirit (though not identical) to the "macroscopic 
compressibility" and "local compressibility" approximations 
introduced (and later abandoned) by these authors (1967A, 

1972) . 

To proceed beyond the lowest order term in f^ - , we note 
again that f 1 is a well-defined, if inaccessible, sum in the 
conventional perturbation series. Perhaps further progress 
can be made by combining our view of f^ with the perturbation 
series approach. As a simple (and possibly important) 
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example, we know from the perturbation series for a pure fluid 

1 

that, whever a single term £v (r) appears, we can perform a 

partial summation of higher order terms by replacing that 

gv^tr) with {l-e”^ v ^}. (See, for example, Brout, 1965, or 

Andersen and Chandler, 1972.) Then we may use (3-75) to gen- 

1 . * 

erate a further approximation to f specified by 


-ijarar' e -6^(£-£') { pP' ref (r ( r') - 

= Jdr dr 1 6p(r)6p(r' ) f 1 (r-r ’ ) + o(6p) 3 . 


( 6 - 4 ) 


(The generalization to multi-component systems is presumably 
straightforward.) The approximation represented by (6-4) must 
represent a valid summation in the perturbation series, and 
goes beyond any we have used in this work. Such an approxima- 
tion will clearly affect the T-dependence of the results. 
(Since, for the metallic systems, we have v 1 explicitly only 
in reciprocal space, the generalization of (3-75) to (6-4) 
introduces non-trivial computational problems.) Direct exten- 
sion of our approach to higher order seems all but out of the 
question, since the expansion leading to (3-75) involves, in 
the next term, the three and four body distribution functions 
of the non - uniform fluid. 

Now consider the adequacy of the overall approach, that 
is, of the approach which calculates the structure by 

*The constant term in {l-e*"^ vl ^ } vanishes from (6-4) by the 
sum rule on p£* (See Appendix A, equation (A-8)). 
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calculating in first order. We note that, in view of the 
objection raised in the discussion of the metallic potentials, 
our metallic calculations may not have provided a test of the 
theory. It may be that the only test we have run on this 
theory is the calculation of S (0) for Na parameters and 
several values of n* While the results of this test were not 
encouraging, it is clearly incomplete and insufficient to 
stand as the only test of the adequacy of the theory. More 
complete and careful testing of this theory on model poten- 
tials is clearly called for. 

Application of this work to Lennard-Jones fluids will 
require a little care, since the work finds its most natural 
expression in reciprocal space. We will thus not be able to 
use the full potential as the perturbation on the hard core as 
we could for metals. Since, however, these results are expli- 
citly independent of the form of v^(r) for r < o , we are at 
liberty to give to v^ a form inside the core which makes 
Fourier transform of v^ possible. This is the approach taken 
by Chandler et al. (1972) in their "optimized" RPA. Direct 
application of this work to molecular liquids must be held 
pending an understanding of the importance of the assumption 
v^fr) = v^(-r), which has apparently been required in several 
places. The essentials of the approach, however, that is, the 
view which takes the densities rather than the potentials as 
the basic variables, may, in suitably generalized form, have 
wide application in the theory of both simple and molecular 
(insulating) liquids. 



-125- 


VI-B. Application to Metallic Liquids 

We have used the results of sections III and IV to define 
a long wavelength structural approximation which we have 
called the mean density approximation (MDA) . In application 
to potentials modelled for the pure alkali metals, this ap- 
proximation yields, at k = 0, results consistent with our 
expectations for the compressibility. At values of k around 
the reference (HC) peak, the results in the MDA seem to give 
some indication of a finite particle size effect. At k = 0 in 
the binary case, a simple model calculation of the line of 
singularities reproduces quite well the asymmetry and position 
of the Li-Na phase boundary in the x-T plane. Because the 
correct line of "hypothetical" singularities lies below the 
phase boundary at all points but the critical point, we expect 
qualitatively the more rounded shape we observe in the calcu- 
lated curve. The selection of the hard sphere reference pa- 
rameters emerges as a central problem. This problem, however, 
was studied almost entirely with respect to metallic poten- 
tials and real systems, so that, in view of the difficulty 
with the metallic potential, we clearly need to do more work 
in the spirit of the test at S{0), that is, in the spirit of 
calculations for model systems. The nature of the singularity 
as k+0 was investigated, and appears asymptotically to be of 
the form (a + bk 2 )~^ (a gratifying result in view of the com- 
plexity of the MDA formulas) , a form which insures the conver- 
gence of the resistivity integral (1-1) at the lower limit. 

The partial structure factors in the MDA were calculated and 
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plotted vs. k for one combination of parameters. Suggestions 
were made for a simple resistivity calculation which might 
shed some light on the nature and origin of the anomaly ob- 
served by Schiirmann and Parks. Finally, a possible difficulty 
of the metallic potential at long wavelength was noted and 
suggestions made for work on the problem. 

These calculations tend to support the hypothesis that 
the separation transition is driven primarily by the tail of 
the potential. This suggests that the study of such phase 
transitions could ultimately yield useful experimental infor- 
mation on this long range part of the potential. Such infor- 
mation has been difficult to extract from the structure fac- 
tors themselves, since these tend to be dominated by the geo- 
metric effects of the repulsive cores. It is obvious, how- 
ever, that an enormous amount of theoretical progress is still 
required. A certain amount of pessimism may be justified. 

A question which remains unanswered is why some metallic 
binaries separate and others do not. One need only consider 
Na-K, rather than Li-Na, to find a binary system in which the 
uniform phase is stable at all concentrations for any temper- 
ature above solidification. Our physical understanding of 
this matter is that phase separation (a long wavelength order- 
ing) and freezing {a short wavelength ordering) must in each 
system be competing transitions. A calculation of sufficient 
subtlety should show that, in Li-Na, the separation transition 
occurs at a higher T than a (hypothetical) freezing transi- 
tion, while the converse is true in Na-K. It remains to be 
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seen whether such a calculation can be performed successfully 
with present techniques. 

The theory of the phase transition implicit in the MDA 
is, of course, essentially mean field, so that we do not ex- 
pect this work to shed any light on the usual problems encoun- 
tered in the study of critical exponents (Fisher, 1964) . The 
critical exponents in this problem are defined by Schurmann 
(1971) , to whose work the interested reader is referred. We 
have paid little attention to this view of the problem and its 
necessary introduction of Ising type models, in favor of an 
approach directly through the structural results and calcula- 
tions of classical liquid theory. 

The structural work has been presented with a resistivity 
calculation ultimately in view. The work, of course, will 
have application near the phase transition to other properties 
of liquid metals which depend on the structure (e.g. the 
thermopower, thermal conductivity, optical absorption, ultra- 
sonic absorption) . Before such problems can be studied with 
any confidence, however, it seems that we must resolve the 
difficulty of the pair potential approximation at long 
wavelength. 



Appendix A: The Distribution Functions 

We shall define the distribution functions for a pure 
fluid. The generalization to mixtures is straightforward but 
cumbersome. We consider a system of N particles in a volume V 
at temperature T, and evaluate thermal averages in a canonical 
ensemble. We define the N particle distribution function 

P n (R> = a - 6U l R > a - M(R >/([dR a a (A-l) 

(just as in (3-9)). This function satisfies 

JdR P n (R) = 1 . (A-2) 

In terms of this function, the thermal average of some quan- 
tity A (R) is calculated, just as in (3-63) , by 

<A(R)> = JdR A(R) P n (R) . (A-3) 

The n-body distribution function is defined by 


‘ > n< r l’ r 2'-" r n ) 


n r r 

TH^HTT J dr n+X--- dr N VW-'-Wn'-'-V • 


(A- 4) 


For n = 1, 


P 1 (r) = N Jd £ 2 • • • dr^ P N (r ,£ 2 , . . . r^) , 


(A- 5) 


is just the single particle density. This function clearly 
satisfies the sum rule 

Jdr P 1 (r) = N . (A-6) 
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For n = 2 , 

@2 ^ = N(N-l) Jdr 3 dr 4 . . .dr^ P N (r,r' ,r 3# . . .r^ ( A " ? ) 

is the two body distribution function which plays the impor- 
tant role in our work. It has the meaning that P 2 (£/£’ )<*£ d £' 
is the probability of finding a particle in dr about r and 
another in dr' about r'. It satisfies the sum rule 


Jdr dr' P 2 (£/£' ) = N(N-l) . 


(A-8) 


Now consider the role of these functions in a thermal 
average, using (A-3) . First consider the thermal average of 
the sum on single particle potentials $(R) defined, as in 
(3-2) , by 


$(R) = JC <*>(£.!> • 


(A-9) 


Then 


<$(R)> » ^ J d £ x d £ 2 * • • d£jj 4>(£ i ) p n (£x , ’**^N ) * ( A_10 > 

Because P N is a symmetric function of its arguments , we may' 
make in each integral in (A-10) the change of variables r^->-r^, 
and r^r^. Si nce all N integrals are then seen to be equiv- 
alent, we have 

<$ (R) > = N^dr^ <f> (£ x ) Jdr 2 . . .d^ P N (r^ ,£ 2 , . . . r^) . (A-ll) 

By (A-5) , this is just 

<<t> (R) > = Jdr ^(rjpj^tr) . (A-12) 
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Next calculate the thermal average of the sum of pair 
potentials U(R) defined, as in (3-42), by 

U(R) = i £ v(r,-r.) . (A-13) 

Then 

<U(R)> = ~ £ | d £i d ^ 2 * ■ * d ^N v( -i“-j^ P K *-l'* * * ( A-14 > 

i^j 

We use the symmetry of P^ (R) to make, in each integral of 
(A-14) , the changes of variable r^->-r^, r^-*-r^ , r^r^. 

Then all N (N-l) integrals in (A-14) are seen to be 
equivalent, and we have 

<U(R)> 

(A-15) 

= | N(N-l) Jdrjdr^ v(r x -r 2 ) P N <£ i /£ 2 . • • •£„) , 

so that, by (A- 7). , we find 

<U(R)> - yjdr dr 1 v(r~r' ) p 2 (£,£' ) • (A-16) 

Finally, consider the thermal average of the density 
operator p (k) defined by (1-3). In real space, this operator 
is 

i3(r) = E«(r-r.) . (A-17) 

i 

Then 

<0(r)> = XZ Jd^dr 2 . . . dr^j 6 (r-r^,) P N (£ 1 ,...r N ) . (A-18) 

Using the symmetry of P^(R), we make, in each integral of 
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( A— 18) , the changes of variable r^-*-r^ and r^+r^, so that all 
N integrals in (A-18) are seen to be equivalent, and we have 


<p (r) > = N J d £x 6 J d £ 2 *** d £N P N ^—1 '—2 ' * ' * 


(A-19) 


= P-l (r) , 


by (A-5) . Then £ (r) is, as stated in section I, the operator 
whose thermal average is the single particle density. Similar 
arguments will establish the relations, given in (1-6) and 
(1-7), between S (k) defined in terms of p (k) (in (1-2)) and 
the two body distribution function p 2 (r ,r ' ) . 



Appendix B: Phase Separation and Singularities in the (k) - 

Here we wish to elaborate on the conclusions stated in 
section I concerning the relation between phase separation and 
singularities in the partial structure factors (k) in the 
limit as k-*0. 

The analysis concerns the form of the Gibbs energy per 
particle, g(x), plotted as a function of x at constant T and 
p. For some < T c , this isotherm will have a form like that 

presented schematically in Figure B-l. This curve exhibits 
regions where it is concave upward and regions where it is 
concave downward. In the region where the curve is concave 
downward, between the points labelled C and D, we have 



so that the system is unstable against concentration fluctua- 
tions. The points A and B are defined (uniquely) by construc- 
ting a line tangent to this curve at two points. A straight- 
forward calculation of the Gibbs energies involved (bearing in 
mind particle number conservation) shows that, for any concen- 
tration between and x g (say x £ ) , the Gibbs energy of the 
uniformly mixed phase is higher than that of some separated 
phase, and that the separated phase with the lowest Gibbs 
energy is that defined by x^ and x fi . On the other hand, for 
concentrations outside this region, x < or x _> Xg, the 
Gibbs energy of the uniform phase is lower than that of any 
separated phase. Then the concentrations and Xg mark the 
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Concentration 


FIGURE B-l 

Gibbs energy per particle as a function of concen- 
tration at T = < T^. 
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points in the x-T plane at which the line T = intersects 
the phase boundary. (Refer to the discussion in section I-A. ) 
As the temperature is raised, these points move closer togeth- 
er, until, at the critical temperature T c , they coincide. For 
all T > T^, the curve is everywhere concave upward. Then a 
common tangent may not be constructed, and the uniform alloy 
is stable in all concentrations. 

Our concern is with the inflection points C and D, at 

which 

= 0 . (B-2 ) 

T,p 

Bhatia and Thornton (1970) have shown that these are the 
points at which the partial structure factors of the uniform 
alloy diverge in the limit as k->-0. It is clear that, at the 
temperature T^ represented by the isotherm in Figure B~l, 

these points are not accessible. They must lie below— the 

phase boundary in the x-T plane. As the temperature is in- 
creased towards T^,, however, these inflection points will move 
closer together until, at T c , the inflection points and the 
common tangent points coincide. It is thus only at the criti- 
cal point that the singularity is fully realized at separation. 
Then the line of singularities in the x-T plane lies below the 
phase boundary in general, but joins it at the one point which 
is the critical point. 




Appendix C: The k = 0 result (4-32). 

We note from (1-2) and (1-3) that 

S..(0) = (N.N.)" 3s <N.N.> - (N.N.)* 5 .. (C-l) 

XJ X J X J X J 

We evaluate <N^N^> in a grand canonical ensemble at constant 
T, V, and chemical potentials vk . Consider the pure fluid 
case, for which 

e^jdR 

<N> . J >— (c-2> 

Ei a^JdR e™ 

N 

Then we find 

(fu) T V <N> = e(<NZ> ~ <N>2 > * (C-3) 

Then 

t ^(w) T jV <n> = ^ <N2> - <N> - s <°> • < c -» 

That is 

s( °> ■ M^) t , v ' where p - • <°- 5 > 

In a multi- component system/ this result generalizes readily 
to 


s 



1 





T,V,y 


/ 

I 


(C-6) 


where the prime indicates that all y's are held fixed except 
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the one explicitly appearing in the derivative. By (3-5) , 
this result is 


Xij (0) 


iM 

^ j/ T , V, y 1 


(C-7) 


Now, if is considered a function of all the y's, and y^ is 
considered a function of all the p's, we have 


-pi - ? 


f»p ± \ 

•r dy . 

i 3y j/T,V,y' 1 


; dy. 


■?&) 


d P, 


(C— 8 ) 


V T,V,p ' 


which gives 



(C-9) 


(C-10) 


or 


P 



(C-ll) 


where P is the matrix whose elements are (3p./3y.) „ ,, and 

1 3 r V / y 

U is the matrix whose elements are (3y./3p 0 )_ . , . Intro- 

ducing the matrices X and F from (3-33), we then have, from 
(C-7) and (C-ll) , 


X = -P = 




(C-12 ) 


so that F = 0, or 
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Frora thermodynamics, this is 


f ij 


(0) = I f 9 2 f - ) 

V l 3 Pi 3 Pj/T,V,p' 


(C-14) 


where F now denotes the free energy of the system. To arrive 
at (4-32) , we expand f^^ and F in powers of the perturbation 
vf j . We introduce an expansion parameter X by replacing v 1 
with Xv 1 , and expand f^j as in section III: 


- f ij f(0> + Xa ij (0 * + o(X)2 


(C-15) 


For the expansion of F, we simply insert the Zwanzig expan 
sion, which, as we saw in the argument to (3-75) , is just 


F = F ref + X<U 1 (R)> ' Tef + o(X) 2 . (C-16) 

X rgf 

On writing <U (R) > - in terms of p?( (r,r')/ just as in 

o r 1 j 

(3-71) , and substituting into (C-14) , the desired result, 
(4-32) , follows immediately. 


We note that this result is also (plausibly) consistent 
with a certain limiting case studied by Lebowitz and Penrose 
(1966). These authors replace v 3 (r) with the form 

v 1 (r) * X 3 <J>(Xr) , (C-17) 

and study the form of the free energy in the thermodynamic 
limit and in the limit as X-*-0. In this limit, the potential 
v 3 becomes very weak and very long range. The limiting result 
for the free energy is 

F(p) = F ref (p) + | P 2 <}>o > 


(C-18) 
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where F{p) now denotes the free energy per unit volume, and 

<f> 0 = Jdr f (r ) =J'dr X 3 <j>(Xr) , for all X. (C-19) 

This result is made plausible by consideration of (C-16) , 
which reads, for this choice of v 1 , 

F(p) = F ref (p) + |Jdr X 3 <^{Xr)p° ' ref (r) + o(cf>) 2 . (C-20) 

(Note that between (C-16) and (C-20) , both F and X have been 
redefined,) If the dimensions of the system are much greater 
than the range of X 3 <t>(Xr), we can change variables in the 
integral to write 

F ( p ) = F ref (p) + |Jdr <|)(r) p° ,ref (r/X) + o(<f>) 2 . (C-21) 

If we now take the thermodynamic limit and the limit as X->-0 
(in that order, so that the dimensions of the system are 
always much greater than the range of X 3 <Ji(Xr)), we see that 
the region of integration in which p^ ' (r/X) differs from 

p 2 will shrink to zero. Then (C-18) will follow if, as seems 
reasonable for this weak form of v"* - (r) , the higher order terms 
may be neglected. 

The result (C-18) yields, through (C-14) , the mean field 
result 

f 1 (0) - <K = Jdr v 1 (r) . (C-22) 

This result can be obtained directly from (4-32) by a similar 
plausibility argument. For this potential, (4-32) reads 
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where, for the last step, we again require that the dimensions 
of the system be much larger than the range of X 3 $(Xr). If we 
perform the same limiting steps described above, the region 
where j p° ,ref (r/X) differs from unity is seen to shrink 

to zero, so that (C-22) follows from the neglect of higher 
order terms. 



Appendix D: The MDA to fT^ (k) for a Two-component Fluid. 

Let x x = p°/p and x 2 = p^/p* where p£ and p° are the den- 
sities of species one and two respectively, and p « p£ + p£. 
Define the partial packing fractions ni and n 2 by , 


where a. is the hard core diameter associated with the i 

i 

species. Then the MDA for f^j for a binary alloy is 


th 


fJl<k) = v^(k) 


(N {n ‘ 


2 157 T - 1 s ? 2 <*-*’> 


+ n.^s^fk-k')) 


^ 4-^ v 

X > 2 PJ (2TT) * 


f 12 v 12 


/ 1, T \ B 1 ^ 

t 22 — ' 2 8rti 2 w 22 v — — 


si!, (k-k * ) , 


(D-l) 


1 1 

r dk' 

x 2 pj 

f (2TT) 


3 2 f k— k 1 } + — 2 — fk— k ' ) } 

Sru3n 2 s ii -- } + 2 3n 2 S ll ( -- )} 


1 if**' 

/XiX 2 P J ( 2 7T ) 


- AL 


3 V 12 (— ' ) ^ nin2 


,HC 


1 „HC 


^ 2 s iJ <*-**> + f 


.L B 1 3 qHC <• u | \ I B2 ft CHC /]_ j- | \ 1 

+ 2T ¥b x S 12 — £ ) + 2 9n 2 S 12 - - U 

v^(k’) |L__ S^(k-k') + §-*- |- S^tk-k')} 

x 1 pj(2 7r ) 3 22 4 — 2 3m3ri2 22 2 3 n x 22 
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The form of f^ 2 (k) may be obtained from the form of f ^ (k) by 
replacing all subscripts ' 1 ' with subscripts '2', and vice 
versa. 

The derivatives of the may be given explicit forms 
in the PY approximation. If these forms are to be derived 
from the published expressions (AL-II) , care must be taken 
that the right variables are held constant as the differen- 
tiations are performed. The structure function S.. is a func- 
tion of k, p°, p°, a and We stipulate, just as in the 

pure fluid case, that the ck are held fixed. Then, for ex- 
ample, some of the necessary derivatives in (4-31) take the 
form 


flfii) 
\ ^1 ' 


(D-2) 


The following discussion of the variables will be valid as 
well for the more complicated derivatives derived from (4-31) , 
but we shall take the expression (D-2) as the simplest example. 

From the definition of the it is clear that we may 

trade the five variables above for k, r^' and a 2 * We 

then have 



Further, since S 
the same amount, 


= n 


,a l ,a 2 


f 8Si , j .) 

1 l 3T1 i /k,n 2 ,a lf ( 


(D-3) 


is unchanged when we scale all lengths by 
we have 


s ij * k ' n i' n 2' a i' a 2* = s i j ^ 1 k/n 1 fn 2 fX<T 1 #Xa 2 ) # (d-4) 


so that may be considered a function of only the four 
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variables ko^, kc? 2 1 and n 2 - Then the expression (D-3) 


becomes 


3S. . 

n 


l\3n 


(D-5) 


1 /ka 1 ,ka 2 ,ri2 


If we now trade ka^ and ka 2 for a = a^/a 2 an< ^ ^°2 ' We arr i ve 


finally at the result 




(D-6) 


where S. . is now considered a function of ka«, a, r\. , and ru. 

ij 2 1 2 

Then the form quoted in AL-I must be rewritten, since the var- 


iables in AL-I are ka 2 , a, x, and n* 
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